# Correct way to populate a DiagonalMatrix?

I would like to create a series of correlation matrices that starts with :

sensMat[[1]] = DiagonalMatrix[ { 1,1,1,1,1 } ]) // MatrixForm


and iterates in 0.1 steps for each of the elements (other than the leading diagonals) until the last matrix of all 1's are obtained. Maybe using the Correlation command?

-
Something like: sensmat = Table[DiagonalMatrix[{1, 1, 1, 1, 1}] /. (0 -> i), {i, 0, 1, 0.1}]; seems simple enough –  ssch Jan 6 '13 at 13:38
Hi @ssh, so the Dimensions[ sensMat[ [ i ;; i] ] ] are {1,5,5}, how to get it to {5,5} as I will use each of the matrices as input argument of {5,5} dimensions for another function. Flatten or Transpose perhaps? –  sebastian c. Jan 6 '13 at 14:42
@sebastianc. why not simply Dimensions[sensMat[[i]]] ? –  P. Fonseca Jan 6 '13 at 15:33
@sebastianc. what @P.Fonseca said, sensmat[[1]] will be the identity matrix and the zeros will increase up until it's all ones in sensmat[[11]]. And @@ (Dimensions[#] == {5, 5} & /@ sensmat) (* True *) –  ssch Jan 6 '13 at 15:57

Five by five matrices can be constructed any way you want--the timing will be fast and space requirements small. But if this is intended to be applied to larger matrices, then the choice of method matters. Built-in matrix operations are usually not only most efficient, but also tend to be expressively clear, as in:

f[r_, n_Integer /; n > 0] := ConstantArray[r, {n, n}] + DiagonalMatrix[ConstantArray[1 - r, n]]


This creates an $n$ by $n$ matrix filled with the constant $r$ and then places $1$s on the diagonal by adding $1-r$ to the diagonal elements only.

The desired sequence of $11$ such matrices is obtained with

f[#, 5] & /@ Range[0, 1, 0.1]


Let's check the timing for larger matrices. This one has $100$ million nonzero entries:

AbsoluteTiming[ArrayPlot[f[0.1, 10^4]];]


{2.6696048, Null}

$60$% of that time is to plot the matrix; the actual computation takes just one second. (The reason for invoking ArrayPlot is that Mathematica is clever about some operations and does not actually perform them until their output is needed: making a plot guarantees the calculation is carried out.)

This is $100$ times faster than, say, iterative solutions or solutions that need to treat the diagonal elements in some special way (conditional expressions--If--are slow).

-

You basically have something analogous to a sparse array with specified elements (the diagonal) in which you want to change the value of the unspecified elements (normally zero). How about this to give you your 11 matrices:

sensMat =Map[SparseArray[{{i_, i_} -> 1}, {5, 5}, #] &, Flatten[{0,Range[0.1, 0.9, 0.1],1}]]


Note that I didn't use Range[0, 1, 0.1] because assume you want 0 and 1 not 0. and 1..

To display:

Normal[sensMat[[3]]]

(* {{1, 0.2, 0.2, 0.2, 0.2}, {0.2, 1, 0.2, 0.2, 0.2}, {0.2, 0.2, 1, 0.2,
0.2}, {0.2, 0.2, 0.2, 1, 0.2}, {0.2, 0.2, 0.2, 0.2, 1}} *)

MatrixForm[sensMat[[3]]]


On my computer this is faster than combining ConstantArray and DiagonalMatrix. The reason this is fast can be found by examining the FullForm of the sparse arrays. You only store one value for unspecified elements. So you store 0 then 0.1 in the next matrix, 0.2 for the next etc. Here is the FullForm of sparse array for the 0.2 matrix above:

SparseArray[Automatic, {5, 5}, 0.2, {1, {{0, 1, 2, 3, 4, 5}, {{1}, {2}, {3}, {4}, {5}}},
{1, 1, 1, 1, 1}}]

-
+1... interesting solution –  P. Fonseca Jan 6 '13 at 21:05
+1 Given that the final matrix is not sparse, it is noteworthy that the SparseArray representation of the diagonal matrix gives a faster calculation. –  whuber Jan 7 '13 at 19:49
@whuber although it is technically not sparse the point was that you can use sparse array functions because all the non diagonal elements are the same value, i.e. they don't have to be zero. –  Mike Honeychurch Jan 7 '13 at 20:49
(1) The SparseArray solution is a lot faster--a couple orders of magnitude. (2) It does not necessarily take much RAM. E.g., making an ArrayPlot with dimensions $10^5$ (involving $10^{10}$ floating point entries) required only about $70$ Mb. (3) If you don't want true zeros or ones at the end of the range, just use Range[0.0, 1.0, 0.1]. If you're concerned about floating point imprecision, Range[0,10]/10.0 should do fine. –  whuber Jan 8 '13 at 13:52
@whuber your first 2 points are the reason for my answer. For (3) I already explained in the answer why I did not use Range[0.0, 1.0, 0.1] -- it was because I mimicked the OPs requirements literally, i.e. he had 0 and 1 (zero and one) in his cut and paste image not 0. and 1.. –  Mike Honeychurch Jan 8 '13 at 19:15

A simple way:

sensMat[n_] := Table[If[i == j, 1, (n - 1)*0.1], {i, 1, 5}, {j, 1, 5}]


A funny way:

sensMat[n_] := NestList[RotateRight, Join[{1}, ConstantArray[(n - 1)*0.1, 4]], 4]


For memorization benefits (and inconvenients):

sensMat[n_] := sensMat[n] = (*your preferred code here*)


-----EDIT-----

Nobody warned me that we were shooting for speed :-), so here is the funny speedy way:

sensMat =
Compile[{{n, _Integer, 0}, {dimension, _Integer, 0}, {step, _Real,
0}}, NestList[RotateRight,
Join[{1}, ConstantArray[(n - 1)*step, dimension - 1]],
dimension - 1]]

AbsoluteTiming[sensMat[3, 10000, 0.1];]
(* {1.670095, Null} *)


In comparison with whuber's solution:

f[r_, n_Integer /; n > 0] :=
ConstantArray[r, {n, n}] + DiagonalMatrix[ConstantArray[1 - r, n]]

AbsoluteTiming[f[0.2, 10000];]
(* {2.828162, Null} *)


(On latter runs, this version gave timings virtually equal to the ones I got with the compiled sensMat)

This last solution doesn't gain anything from compilation (as expected...):

fCom = Compile[{{r, _Real, 0}, {n, _Integer, 0}},
ConstantArray[r, {n, n}] +
DiagonalMatrix[ConstantArray[1 - r, n]]];

AbsoluteTiming[fCom[0.2, 10000];]
(* {2.840162, Null} *)


Just to make sure:

sensMat[3, 10000, 0.1] == fCom[0.2, 10000]
(* True *)


Interesting enough, this is the first case I found where WVM is faster than C

-
Wouldn't you want to use memorization so that matrices don't have to be constructed each time they are called? –  David Carraher Jan 6 '13 at 16:34
@dude added it. Thanks –  P. Fonseca Jan 6 '13 at 16:46
+1 for the edit. But have you honestly compared the timing by forcing an evaluation of both results, such as with an ArrayPlot? When I do that, I find the compiled solution is only 25% faster. That's a nice achievement, but because it is so little, I would in many cases prefer a clear simple native Mathematica solution to a compiled (and perhaps obscure) solution. As you hinted, speed isn't everything! –  whuber Jan 7 '13 at 20:01
@whuber I agree that your solution is more "balanced" –  P. Fonseca Jan 8 '13 at 11:57