Re: Kronecker Product with SparseArray
- To: mathgroup at smc.vnet.net
- Subject: [mg62144] Re: [mg62108] Kronecker Product with SparseArray
- From: "Carl K. Woll" <carl at woll2woll.com>
- Date: Sat, 12 Nov 2005 03:32:44 -0500 (EST)
- References: <200511110752.CAA29703@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Cyrus F Hirjibehedin wrote:
> To Mathgroup:
>
> I'm trying to implement a function in Mathematica 5.2 to take the
> Kronecker Product of two large sparse matrices. The simplest suggested
> way to do this seems to be
>
> BlockMatrix[Outer[Times, X, Y]]
>
> where X and Y are two matrices. The BlockMatrix is necessary because
> Outer produces a matrix of matrices, rather than a matrix of elements, and
> this is not compatible with functions like Eigenvalues or Eigensystem.
>
> What I would like to do is pass in X and Y as SparseArray types and have
> the result produce a SparseArray type that I can pass to functions like
> Eigenvalues or Eigensystem. Right now, I'm doing this with
>
> SparseArray[BlockMatrix[Outer[Times, X, Y]]]
>
> where X and Y are SparseArray types.
An alternative is to use a blockmatrix function which doesn't break up
sparse arrays (so that the extra SparseArray call is unnecessary):
bm[b_] := Join @@ (Transpose[Join @@ Transpose /@ #] &) /@ b
Unfortunately this doesn't seem to help your memory problem.
> The problem is that I run out of
> memory when trying to compute large matrices. I think the problem is that
> when BlockMatrix tries to "flatten" the results of Outer it first
> generates the full dense matrix representation, which is too big either
> for it or for SparseArray to handle.
>
> Does anyone know how to turn the results of Outer[Times,X,Y]] into an
> object that is a SparseArray type that can be passed to things like
> Eigenvalues or Eigensystem?
>
> or
>
> Does anyone know of a different way to directly compute the Kronecker
> Product of two matrices represented as SparseArray types and have the
> output be a SparseArray type that can be passed to Eigenvalues without
> generating an intermediate dense matrix representation?
>
One idea is to avoid using Outer, which constructs a rank 4 tensor which
you have to then convert to a rank 2 tensor. Instead, extract the
nonzero elements and their positions from the sparse arrays, use this
information to create a list of the nonzero elements and positions of
the kronecker product, and then construct the sparse array. The
following function does this:
Off[General::spell1];
Off[General::spell];
toPA = Developer`ToPackedArray;
kroneckerproduct[m_, n_] := Module[{mrules, nrules, mind1, mind2, mval,
nind1, nind2, nval, mnind, mnval},
mrules = Most[ArrayRules[m]];
nrules = Most[ArrayRules[n]];
mind1 = Dimensions[n][[1]](toPA[mrules[[All, 1, 1]]] - 1);
mind2 = Dimensions[n][[2]](toPA[mrules[[All, 1, 2]]] - 1);
mval = toPA[mrules[[All, 2]]];
nind1 = toPA[nrules[[All, 1, 1]]];
nind2 = toPA[nrules[[All, 1, 2]]];
nval = toPA[nrules[[All, 2]]];
mnind = Transpose[{
Flatten[toPA[Outer[Plus, mind1, nind1]]],
Flatten[toPA[Outer[Plus, mind2, nind2]]]
}];
mnval = Flatten[toPA[Outer[Times, mval, nval]]];
SparseArray[mnind -> mnval, Dimensions[m]Dimensions[n]]
]
See the mathworld entry on matrix direct products for an explanation of
the mind1 and mind2 lines and some of the algebra. In my tests with 5.2,
this function seems to use much less memory. For example:
SeedRandom[1];
{a, b} = SparseArray /@ Table[Random[Integer], {2}, {100}, {50}];
With this data, bm[Outer[Times,a,b]] had a MaxMemoryUsed[] of
~325Mbytes, while kroneckerproduct[a,b] had a MaxMemoryUsed[] of
~208Mbytes. Note that just Outer[Times,a,b] had a MaxMemoryUsed[] of
~257Mbytes.
> Any advice would be greatly appreciated.
>
> Sincerely,
>
> Cyrus
>
> --
>
> Cyrus F. Hirjibehedin
> hirjibe at us.ibm.com
Carl Woll
Wolfram Research
- References:
- Kronecker Product with SparseArray
- From: Cyrus F Hirjibehedin <hirjibe@us.ibm.com>
- Kronecker Product with SparseArray