Alex Rivera | Logout

What is the most efficient way to construct large block matrices in Mathematica?

Asked 2011-07-28T23:40:42.963
11

Inspired by Mike Bantegui's question on constructing a matrix defined as a recurrence relation, I wonder if there is any general guidance that could be given on setting up large block matrices in the least computation time. In my experience, constructing the blocks and then putting them together can be quite inefficient (thus my answer was actually slower than Mike's original code). Join and possibly ArrayFlatten are possibly less efficient than they could be.

Obviously if the matrix is sparse, one can use SparseMatrix constructs, but there will be times when the block matrix you are constructing is not sparse.

What is best practice for this kind of problem? I am assuming the elements of the matrix are numeric.

matrix wolfram-mathematica
Edit
Report

1 Answer

8

The code shown below is available here: http://pastebin.com/4PWWxGhB. Just copy and paste it into a notebook to test it out.

I was actually trying to do several functional ways of calculating matrices, since I figured the functional way (which is typically idiomatic in Mathematica) is more efficient.

As one example, I had this matrix which was composed of two lists:

In: L = 1200;
e = Table[..., {2L}];
f = Table[..., {2L}];

h = Table[0, {2L}, {2L}];
Do[h[[i, i]] = e[[i]], {i, 1, L}];
Do[h[[i, i]] = e[[i-L]], {i, L+1, 2L}];
Do[h[[i, j]] = f[[i]]f[[j-L]], {i, 1, L}, {j, L+1, 2L}];
Do[h[[i, j]] = h[[j, i]], {i, 1, 2 L}, {j, 1, i}];

My first step was to time everything.

In: h = Table[0, {2 L}, {2 L}];
AbsoluteTiming[Do[h[[i, i]] = e[[i]], {i, 1, L}];]
AbsoluteTiming[Do[h[[i, i]] = e[[i - L]], {i, L + 1, 2 L}];]
AbsoluteTiming[
 Do[h[[i, j]] = f[[i]] f[[j - L]], {i, 1, L}, {j, L + 1, 2 L}];]
AbsoluteTiming[Do[h[[i, j]] = h[[j, i]], {i, 1, 2 L}, {j, 1, i}];]

Out: {0.0020001, Null}
{0.0030002, Null}
{5.0012861, Null}
{4.0622324, Null}

DiagonalMatrix[...] was slower than the do loops, so I decided to just use Do loops on the last step. As you can see, using Outer[Times, f, f] was much faster in this case.

I then wrote the equivalent using Outer for the blocks in the upper right and bottom left of the matrix, and DiagonalMatrix for the diagonal:

AbsoluteTiming[h1 = ArrayPad[Outer[Times, f, f], {{0, L}, {L, 0}}];]
AbsoluteTiming[h1 += Transpose[h1];]
AbsoluteTiming[h1 += DiagonalMatrix[Join[e, e]];]


Out: {0.9960570, Null}
{0.3770216, Null}
{0.0160009, Null}

The DiagonalMatrix was actually slower. I could replace this with just the Do loops, but I kept it because it was cleaner look

answered 2011-07-29T00:54:08.660

Your Answer