Re: Kronecker Product with SparseArray
- To: mathgroup at smc.vnet.net
- Subject: [mg62147] Re: Kronecker Product with SparseArray
- From: Peter Pein <petsie at dordos.net>
- Date: Sat, 12 Nov 2005 03:32:53 -0500 (EST)
- References: <dl1k06$32$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Cyrus F Hirjibehedin schrieb:
> 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. 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?
>
> Any advice would be greatly appreciated.
>
> Sincerely,
>
> Cyrus
>
> --
>
> Cyrus F. Hirjibehedin
> hirjibe at us.ibm.com
Hello Cyrus,
one approach is to manipulate the ArrayRules directly:
In[1]:=
directProduct[m1_SparseArray, m2_SparseArray] :=
Module[{r, c, pru},
{r, c} = Dimensions[m2];
pru = ArrayRules[m1] /. HoldPattern[{i_Integer, j_Integer} -> v_] :>
Cases[ArrayRules[m2], HoldPattern[{k_Integer, l_Integer} -> w_] :>
{r, c}*{i - 1, j - 1} + {k, l} -> v*w];
SparseArray[Flatten[pru], Dimensions[m1]*{r, c}]]
In[2]:=
A = SparseArray[{{i_, j_} -> a[i, j]}, {2, 2}];
B = SparseArray[{{i_, j_} -> b[i, j]}, {3, 2}];
In[4]:=
directProduct[A, B]
Out[4]=
SparseArray[<24>,{6,4}]
In[5]:=
Normal[%]
Out[5]=
{{a[1, 1]*b[1, 1], a[1, 1]*b[1, 2], a[1, 2]*b[1, 1], a[1, 2]*b[1, 2]},
{a[1, 1]*b[2, 1], a[1, 1]*b[2, 2], a[1, 2]*b[2, 1], a[1, 2]*b[2, 2]},
{a[1, 1]*b[3, 1], a[1, 1]*b[3, 2], a[1, 2]*b[3, 1], a[1, 2]*b[3, 2]},
{a[2, 1]*b[1, 1], a[2, 1]*b[1, 2], a[2, 2]*b[1, 1], a[2, 2]*b[1, 2]},
{a[2, 1]*b[2, 1], a[2, 1]*b[2, 2], a[2, 2]*b[2, 1], a[2, 2]*b[2, 2]},
{a[2, 1]*b[3, 1], a[2, 1]*b[3, 2], a[2, 2]*b[3, 1], a[2, 2]*b[3, 2]}}
In[6]:=
SeedRandom[1];
bigA =
SparseArray[{i_, j_} /; Abs[i - j] == 1 :> 2*Random[] - 1, {300, 300}]
Out[7]=
SparseArray[<598>,{300,300}]
In[8]:=
Timing[mp = directProduct[bigA, bigA]]
Out[8]=
{3.734*Second, SparseArray[<357604>,{90000,90000}]}
In[9]:=
Timing[Chop[Eigenvalues[mp, 10]]]
Out[9]=
{14.672*Second,
{-1.3620353730438348, 1.3620353730435342, -1.3620353730435304,
1.3620353730435268, 1.3207994704176738*I, -1.3207994704176738*I,
-1.3060063042433008, 1.306006304243299, 1.3050519740014384*I,
-1.3050519740014384*I}}
This function directProduct[] works only for SparseArrays with default
value 0.
Peter