3
$\begingroup$

Is there a fast way to construct this matrix?

$\left[\begin{array}{c} I_4\otimes e_1\\ \vdots \\ I_4 \otimes e_T \end{array}\right]$

$e_i$ is the $i$-th column of the matrix $I_T$, $\otimes$ is the Kronecker product.

$I_4\otimes e_t= \left[\begin{array}{cccc} e_t & 0\cdot e_t & 0 \cdot e_t & 0 \cdot e_t\\ 0 \cdot e_t & e_t & 0 \cdot e_t & 0 \cdot e_t\\ 0 \cdot e_t & 0 \cdot e_t & e_t & 0 \cdot e_t\\ 0 \cdot e_t & 0 \cdot e_t & 0 \cdot e_t & e_t \end{array}\right]_{4T\times 4}$

$e_t$ has dimensions $T\times 1$.

$T$ will be in the order of at least 200.

I'm looking for a "vectorized way" of building this matrix. I don't want to have to define it elementwisely, since this will take a lot of time for big values of $ T $.

$\endgroup$
  • 2
    $\begingroup$ Have you tried KroneckerProduct? $\endgroup$ – bill s Nov 13 '18 at 14:42
  • $\begingroup$ @bills I've edited my question. I know the existence of that function. My problem is that I want to create the matrix as a whole, without having to define it componentwise... $\endgroup$ – An old man in the sea. Nov 13 '18 at 14:54
  • $\begingroup$ So what Dimensions do the inputs get and what Dimensions should the output have? $\endgroup$ – Αλέξανδρος Ζεγγ Nov 13 '18 at 16:05
  • $\begingroup$ @ΑλέξανδροςΖεγγ I've added the dimensions. The only thing that may change is the value of $T$ depending on the iteration. $\endgroup$ – An old man in the sea. Nov 13 '18 at 22:33
  • $\begingroup$ Is $ I_T $ the $ T $ dimensional identity matrix? $\endgroup$ – Αλέξανδρος Ζεγγ Nov 14 '18 at 2:41
1
$\begingroup$

Update You may use the ArrayRules syntax of SparseArray.

ClearAll[specialFromRules]
specialFromRules[t_, i_] :=
 SparseArray[
  Flatten@
   MapIndexed[
    {(#1 - 1) (t i) + #2 t - (t - #1), #2} & @@ #2 -> #1 &,
    ConstantArray[1, {t, i}],
    {2}
    ],
  {t t i , i}
  ]

then

specialFromRules[2, 2] // MatrixForm

Mathematica graphics

and

specialFromRules[2, 4] // MatrixForm

Mathematica graphics



Previous method from IdentityMatrix is not as fast as above SparseArray ArrayRules method.

You may use SparseArray with Band.

With a helper function to build the blocks,

ClearAll[specialMatrix]
specialMatrix[t_, i_] :=
 Function[{et},
   SparseArray[
    Map[
     Band[{t (# - 1) + 1 , #}, {t #, #}, {1, 0}] -> et &,
     Range@i
     ],
    {i*t, i}
    ]
   ] /@ Transpose@IdentityMatrix[t]

Then

ArrayFlatten[ List /@ specialMatrix[2, 2]] // MatrixForm

Mathematica graphics

and

ArrayFlatten[ List /@ specialMatrix[2, 4]] // MatrixForm

Mathematica graphics

and so on.

Hope this helps.

$\endgroup$
  • $\begingroup$ Edmund, thank you for your answer and +1. I was wondering if sparse matrix representation is the most efficient way to do computations with the big matrix. Also, can the sparse matrix in mathematica be used to represent $I_{T}\otimes K$, where K has dim $T \times T$ and all entries non-null? $\endgroup$ – An old man in the sea. Nov 14 '18 at 13:33
2
$\begingroup$

I guess that the obvious (as far as I'm concerned) way to go about solving the problem at hand looks something like that:

With[{i4 = IdentityMatrix[4]},

  f[T_] := With[{iT = IdentityMatrix[T]},
    Flatten[KroneckerProduct[i4, #] & /@ iT, 1]
   ]

 ]

Taking into consideration the special structure of the desired matrix $\left[\begin{array}{c} I_4\otimes e_1\\ \vdots \\ I_4 \otimes e_T \end{array}\right]$ we can define another function customKroneckerProduct:

With[{a = {{1, 1}, {2, 3}, {3, 5}, {4, 7}}, b = {{0, 0}, {0, 1}, {0, 2}, {0, 3}}},

  Clear[init];

  init[T_] := a + b (T - 2)

 ];

Block[{plastic},

  ClearAll[plus];
  SetAttributes[plus, Listable];

  plus[plastic[x__]] := List[x] + {4, 1};

  setup[T_] := NestList[
    plus[Apply[plastic, #, 1]] &, init[T], T - 1];

  customKroneckerProduct[T_] := SparseArray[Join @@ setup[T] -> 1, {T, T} 4]
 ]

It seems that customKroneckerProduct returns in approximately half the time-or even faster-it takes f to produce its results, predominantly for problem instances of large size:

enter image description here

$\endgroup$
  • $\begingroup$ Thanks for the answer. Unfortunately, the answer is not correct I think... I've added some more info on the question. Thank you for being willing to help me. $\endgroup$ – An old man in the sea. Nov 13 '18 at 22:30
  • $\begingroup$ @Anoldmaninthesea. you're really kind. I read mistakenly the subscripts in your question but I am confident the answer can be adapted; will try to as soon as possible; thanks! $\endgroup$ – user42582 Nov 14 '18 at 7:17
  • $\begingroup$ +1 and thanks for the answer. ;) $\endgroup$ – An old man in the sea. Nov 14 '18 at 13:29
  • $\begingroup$ you're welcome, thanks for the nice question $\endgroup$ – user42582 Nov 14 '18 at 13:38
2
$\begingroup$

Updated

(the function construct2 has been updated to be simpler and a bit faster)

The natural method to construct your matrix:

construct[i4_, it_] := Join @@ Table[
    KroneckerProduct[i4, Take[it, All, {k}]],
    {k, Length[it]}
]

Here's an alternative that does a single KroneckerProduct, and then massages it's shape:

construct2[i4_,it_]:=With[{len=Length[i4],t=Length[it]},
    ArrayReshape[
        Transpose @ ArrayReshape[
            KroneckerProduct[i4, Transpose @ it],
            {len, len t^2}
        ],
        {len t^2,len}
    ]
]

Simple example:

r1 = construct[IdentityMatrix[4], Array[Subscript[a, ##]&, {2,2}]];
r2 = construct2[IdentityMatrix[4], Array[Subscript[a, ##]&, {2,2}]];

r1 === r2

r1 //MatrixForm //TeXForm

True

$\left( \begin{array}{cccc} a_{1,1} & 0 & 0 & 0 \\ a_{2,1} & 0 & 0 & 0 \\ 0 & a_{1,1} & 0 & 0 \\ 0 & a_{2,1} & 0 & 0 \\ 0 & 0 & a_{1,1} & 0 \\ 0 & 0 & a_{2,1} & 0 \\ 0 & 0 & 0 & a_{1,1} \\ 0 & 0 & 0 & a_{2,1} \\ a_{1,2} & 0 & 0 & 0 \\ a_{2,2} & 0 & 0 & 0 \\ 0 & a_{1,2} & 0 & 0 \\ 0 & a_{2,2} & 0 & 0 \\ 0 & 0 & a_{1,2} & 0 \\ 0 & 0 & a_{2,2} & 0 \\ 0 & 0 & 0 & a_{1,2} \\ 0 & 0 & 0 & a_{2,2} \\ \end{array} \right)$

Timing comparison:

m = RandomReal[1, {2000, 2000}];

r1 = construct[IdentityMatrix[4], m]; //AbsoluteTiming
r2 = construct2[IdentityMatrix[4], m]; //AbsoluteTiming

r1 === r2

{5.97644, Null}

{0.725066, Null}

True

Addendum

The OP asked about returning a SparseArray representation instead. Since most Mathematica functions work seamlessly with SparseArray objects, you just need to use SparseArray objects as arguments to construct or construct2 to get out a SparseArray object. Consider the following two example $I_T$ matrices:

dense = RandomReal[1, {2000, 2000}];
sparse = SparseArray[Band[{1,1}]->RandomReal[1, 2000]];

Then for the dense matrix second argument:

d1 = construct[IdentityMatrix[4, SparseArray], dense]; //AbsoluteTiming
d2 = construct2[IdentityMatrix[4, SparseArray], dense]; //AbsoluteTiming
Head /@ {d1, d2}
d1 === d2

{0.906973, Null}

{1.42632, Null}

{SparseArray, SparseArray}

True

And for the sparse second matrix second argument:

d1 = construct[IdentityMatrix[4, SparseArray], sparse]; //AbsoluteTiming
d2 = construct2[IdentityMatrix[4, SparseArray], sparse]; //AbsoluteTiming
Head /@ {d1, d2}
d1 === d2

{0.314939, Null}

{0.087633, Null}

{SparseArray, SparseArray}

True

Summarizing, if both matrices are sparse, then use construct2, but if the second matrix is dense, use construct. In both cases the output will be a sparse array.

$\endgroup$
  • $\begingroup$ Carl, thank you for your answer. I like it. However, I'm not sure I should implement it, instead of sparse matrices, since I'll have to do some computations with the this big matrix. What do you recommend? $\endgroup$ – An old man in the sea. Nov 14 '18 at 13:28
  • $\begingroup$ @Anoldmaninthesea. Please see update. $\endgroup$ – Carl Woll Nov 14 '18 at 15:49
  • $\begingroup$ +1 thanks Carl ;) $\endgroup$ – An old man in the sea. Nov 15 '18 at 13:02
0
$\begingroup$

I must admit, I am not entirely sure what you are looking for, but might this be along the right lines?

I4 = {{i11, i12, i13, i14}, {i21, i22, i23, i24}, {i31, i32, i33, 
i34}, {i41, i42, i43, i44}}; MatrixForm@I4

I4 matrix

Transpose[KroneckerProduct[I4, #] & /@ Transpose[I4]] // MatrixForm

enter image description here

$\endgroup$
  • $\begingroup$ Thanks for the answer. Unfortunately, the answer is not correct I think... I've added some more info on the question. Thank you for the help. $\endgroup$ – An old man in the sea. Nov 13 '18 at 22:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.