Re: DirectSum (feature request)
- To: mathgroup at smc.vnet.net
- Subject: [mg97818] Re: DirectSum (feature request)
- From: Raffy <raffy at mac.com>
- Date: Sun, 22 Mar 2009 05:50:17 -0500 (EST)
- References: <gq2eu9$ech$1@smc.vnet.net>
On Mar 21, 3:18 am, Maris Ozols <maroz... at gmail.com> wrote: > Taking a direct sum of a given list of matrices is a very common task > (unless you are a quantum physicist and use only KroneckerProduct). > Unfortunately there is no built-in function (that I know of) for doing > this in Mathematica. The closest thing we have is ArrayFlatten. So I > usually do something like this to compute a direct sum: > > DirectSum[Ms_] := Module[{n = Length[Ms], z, i}, > z = ConstantArray[0, n]; > ArrayFlatten@Table[ReplacePart[z, i -> Ms[[i]]], {i, 1, n}] > ]; > > Is there a better way of doing this? > > Note: A nice way to implement it would be > > DirectSum[Ms_] := ArrayFlatten@DiagonalMatrix[Ms]; > > Unfortunately this gives "DiagonalMatrix::vector" error, since > DiagonalMatrix is not flexible enough to accept a list of matrices. > The way DiagonalMatrix is used in the above code might cause some > confusion for beginners, but in general I don't see why DiagonalMatrix > should be limited in this way. > > ~Maris Ozols~ overhead[_?(VectorQ[#, MatrixQ] &)] = 0; directSum1[vComp_?(VectorQ[#, MatrixQ] &)] := Block[ {mForm = ConstantArray[0, {Length[vComp], Length[vComp]}]}, Do[mForm[[i, i]] = vComp[[i]], {i, Length[vComp]}]; ArrayFlatten[mForm]]; directSum2[vComp_?(VectorQ[#, MatrixQ] &)] := ArrayFlatten@ Array[If[Equal[##], vComp[[#1]], 0] &, {Length[vComp], Length[vComp]}]; directSum3[vComp_?(VectorQ[#, MatrixQ] &)] := ArrayFlatten[ DiagonalMatrix[Range@Length[vComp]] /. i_ /; i > 0 :> vComp[[i]]]; directSum4[vComp_?(VectorQ[#, MatrixQ] &)] := ArrayFlatten@ReleaseHold@DiagonalMatrix[Hold /@ vSum]; directSum5[vComp_?(VectorQ[#, MatrixQ] &)] := ArrayFlatten[DiagonalMatrix[$ /@ vSum] /. $ -> Sequence]; vSum= {{{a, b}, {c, d}}, {{n}}, {{r, s, t}, {u, v, w}, {x, y, z}}}; In[1718]:= Do[overhead[vSum], {10000}] // Timing Out[1718]= {0.025205, Null} In[1721]:= Do[r1 = directSum1[vSum], {10000}] // Timing Out[1721]= {0.751228, Null} In[1727]:= Do[r2 = directSum2[vSum], {10000}] // Timing Out[1727]= {0.79516, Null} In[1731]:= Do[r3 = directSum3[vSum], {10000}] // Timing Out[1731]= {0.857756, Null} In[1749]:= Do[r4 = directSum4[vSum], {10000}] // Timing Out[1749]= {0.65793, Null} In[1764]:= Do[r5 = directSum5[vSum], {10000}] // Timing Out[1764]= {0.71162, Null} directSum4 seems to be the fastest, it is also the simplest and most straightforward. The Hold/ReleaseHold is used to cheat the DiagonalMatrix VectorQ check.