I have a matrix that I generate with a procedure that looks like this:
blaMX[ptsX_Integer, ptsY_Integer] :=
Table[
With[
{
ix = 1 + Floor[(i - 1)/(ptsY)], jx = 1 + Floor[(j - 1)/(ptsY)],
iy = Mod[i, ptsY, 1], jy = Mod[j, ptsY, 1]
},
If[iy != jy,
0,
k1[ix, jx]
] +
If[ix != jx,
0,
k2[iy, jy]
]
],
{i, ptsX*ptsY},
{j, ptsX*ptsY}
];
where I have two square matrices k1 of dimension ptsX and k2 of dimension ptsY
It generates a block-diagonal matrix with a bunch of block-bands. The basic structure looks like:
blaMX[6, 9] /. {_k1 -> 1, _k2 -> 2, c -> 3} // MatrixPlot
where each of those diagonal blocks is k2 duplicated:
blaMX[3, 6][[7 ;; 12, 7 ;; 12]] /. _k1 -> 0 // Grid
k2[1,1] k2[1,2] k2[1,3] k2[1,4] k2[1,5] k2[1,6]
k2[2,1] k2[2,2] k2[2,3] k2[2,4] k2[2,5] k2[2,6]
k2[3,1] k2[3,2] k2[3,3] k2[3,4] k2[3,5] k2[3,6]
k2[4,1] k2[4,2] k2[4,3] k2[4,4] k2[4,5] k2[4,6]
k2[5,1] k2[5,2] k2[5,3] k2[5,4] k2[5,5] k2[5,6]
k2[6,1] k2[6,2] k2[6,3] k2[6,4] k2[6,5] k2[6,6]
Each band is made up of blocks, as an unwrapping of k1 with the Band starting at {n+(i-1)*ptsY, n+(j-1)*ptsY} being k1[[i+1, j+1]] (essentially just using Quotient:
With[{n = 1},
Table[
blaMX[6, 9][[
(i - 1)*9 + n,
(j - 1)*9 + n
]] /. _k2 -> 0,
{i, 6},
{j, 6}
]
] // Grid
Now, my generation procedure is wildly inefficient. I do ptsX^2 * ptsY^2 work, when I really only have ptsX * ptsY * (ptsX + ptsY - 1) (which is just 2 n^3 - n^2 if ptsX == ptsY == n) non-zero elements.
So, what is the best way to build a matrix of this form? My current version looks like:
makeKMat[k1_, k2_] :=
With[{
k1SparseRules =
With[{ptsX = Length@k1, ptsY = Length@k2},
Flatten[
Table[
Band[{(i - 1)*ptsY + 1, (j - 1)*ptsY + 1}, {i*ptsY, j*ptsY}] ->
k1[[i, j]],
{i, ptsX},
{j, ptsX}
],
1
]
],
k2SparseRules =
With[{ptsX = Length@k1, ptsY = Length@k2},
{
Band[{1, 1}] ->
ConstantArray[k2, ptsX]
}
]
},
With[{k3 = SparseArray[k2SparseRules]},
SparseArray[k1SparseRules, Dimensions[k3]] + k3
]
];
And this works:
k1Pts = 6;
k2Pts = 9;
k1Test = Array[k1, {k1Pts, k1Pts}];
k2Test = Array[k2, {k2Pts, k2Pts}];
kMatTest = makeKMat[k1Test, k2Test];
kMatTest == blaMX[6, 9]
True
And is pretty fast, all considered:
k1Test2 = RandomReal[{}, {15, 15}];
k2Test2 = RandomReal[{}, {15, 15}];
{tt, mx} = makeKMat[k1Test2, k2Test2] // RepeatedTiming;
tt
0.0071
{tt2, lmx} = lazyMX[k1Test2, k2Test2] // RepeatedTiming;
tt2
0.338
But, I can't help but feel that this could be done more efficiently. Can someone propose a better way to build a block-banded SparseArray like this?
Some mild goal-post moving:
What if we wanted to extend this to an $N$-dimensional system? A way to implement this element-wise would be:
getCrds[i_, j_, gridDivs_] :=
FoldList[
With[{tot = #[[3]], new = #2},
Append[
Map[
Mod[1 + Quotient[# - 1, tot], new, 1] &,
{i, j}
],
tot*new
]
] &,
{i, j, 1},
gridDivs
][[2 ;;, ;; 2]];
makeNDKMat[ks : {__List}] :=
With[{gridDivs = Length /@ ks},
Table[
With[{divvy = getCrds[i, j, gridDivs]},
Total@
Table[
If[AllTrue[Delete[divvy, k], Apply[Equal]],
ks[[k, Sequence @@ divvy[[k]]]],
0
],
{k, Length@gridDivs}
]
],
{i, Times @@ gridDivs},
{j, Times @@ gridDivs}
]
]
For example, in 4D this looks like:
k1T = Array[k1, {2, 2}];
k2T = Array[k2, {3, 3}];
k3T = Array[k3, {4, 4}];
k4T = Array[k4, {5, 5}];
k4DLazy =
makeNDKMat[{k4T, k3T, k2T, k1T}] /. {_k1 -> 1, _k2 -> 2, _k3 ->
3, _k4 -> 4};
k4DLazy // MatrixPlot
How could we efficiently build a sparse representation of this?
