MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Kronecker Product with SparseArray

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:

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:

{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 

> Any advice would be greatly appreciated.
> Sincerely,
> Cyrus
> --
> Cyrus F. Hirjibehedin
> hirjibe at

Carl Woll
Wolfram Research

  • Prev by Date: Re: How long does technical support take?
  • Next by Date: Re: How long does technical support take?
  • Previous by thread: Kronecker Product with SparseArray
  • Next by thread: Re: Kronecker Product with SparseArray