Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2005
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2005

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

Search the Archive

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


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