Fast way to create $[ I_4\otimes e_1,\ \dots ,\ I_4 \otimes e_T]$?

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$$.

• Have you tried KroneckerProduct? Nov 13 '18 at 14:42
• @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... Nov 13 '18 at 14:54
• So what Dimensions do the inputs get and what Dimensions should the output have? Nov 13 '18 at 16:05
• @ΑλέξανδροςΖεγγ I've added the dimensions. The only thing that may change is the value of $T$ depending on the iteration. Nov 13 '18 at 22:33
• Is $I_T$ the $T$ dimensional identity matrix? Nov 14 '18 at 2:41

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


and

specialFromRules[2, 4] // MatrixForm


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


and

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


and so on.

Hope this helps.

• 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? Nov 14 '18 at 13:33

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:

• 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. Nov 13 '18 at 22:30
• @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! Nov 14 '18 at 7:17
• +1 and thanks for the answer. ;) Nov 14 '18 at 13:29
• you're welcome, thanks for the nice question Nov 14 '18 at 13:38

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

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
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
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.

• 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? Nov 14 '18 at 13:28
• @Anoldmaninthesea. Please see update. Nov 14 '18 at 15:49
• +1 thanks Carl ;) Nov 15 '18 at 13:02

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


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


• 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. Nov 13 '18 at 22:30