How can I define a very large matrix efficiently?

size = 7;
Table[Which[i == 1, 1, i > j, 1, i == j, -i + 1, True, 0], {i, size}, {j, size}]

The question refers to efficiency, which might mean with respect to programming, performance, or both. My guess is programming, but highlighting the performance strengths of Mathematica should promote learning the "Mathematica style" as it is sometimes referred to.

Timings on various approaches for size = 100 (timing code given at the end):

  • 0.008800 @kglr
  • 0.007000 @Rohit Namjoshi, which can be improved to 0.00026
  • 0.000510 @MikeY, which can be improved to 0.000041 (best of all)
  • 0.000047 @MichaelE2, (was never under 0.000042 but usually 0.000046 to 0.000048)

My best solution:

size = 5;
matrix = 
 With[{r    = Range@size,
       zero = Developer`ToPackedArray@{0}},
    UnitStep[Outer[Plus, r, -r] + PadRight[ConstantArray[r, {1}], {size, size}]] - 
     DiagonalMatrix@PadRight[zero, size, r]
    ];
matrix // MatrixForm

Mathematica graphics

Discussion

Table[] vs. Array[]. The solution of @RohitNamjoshi can be improved by using Array[]. This involves turning the expression e into a function, which takes only wrapping the expression with Function[{i, j}, e] and removing i and j from the iterators:

Array[Function[{i, j},
    Which[i == 1, 1, i > j, 1, i == j, -i + 1, True, 0]
    ], {size, size}]; // RepeatedTiming
(*  {0.00039, Null}  *)

It can be sped up a little by compiling the function. Compiling with or without CompilationTarget -> "C" affects Table[] and Array[] differently. Compiling speeds up Table[] (0.0039 sec.), but whether to C or not makes no perceptible difference. Compiling to C speeds up Array[] slightly, but compiling to WVM the option slows down Array[] (0.00091 sec.); I don't know why.

cf = Compile[{{i, _Integer}, {j, _Integer}},
   Which[
    i == 1, 1,
    i > j, 1,
    i == j, -i + 1,
    True, 0], CompilationTarget -> "C"];

Array[cf, {size, size}]; // RepeatedTiming
(*  {0.00026, Null}  *)

Vectorization vs. Compiling. Many basic functions in Mathematica are vectorized: They work efficiently on packed arrays. Elementary functions, many linear algebra functions, some list manipulation functions (esp. functions operating on rectangular arrays) are particularly efficient on packed arrays. Other functions might "unpack" the array, which involves copying the packed array into the unpacked form. This done internally but there are ways to tell (On["Packing"] and Developer`PackedArrayQ for instance); mainly it is invisible to the user except for a performance hit. Usually, but not always, the vectorized functions are more efficient than compiled functions. The solution of @MikeY takes advantage of this up until ones and the use of ArrayFlatten[], which unpacks its arguments to assemble the final matrix. The reason ArrayFlatten[] unpacks is because one of the inputs, ones, is not packed. It turns out that user-entered lists are not packed. If we pack it, we get a great improvement in performance:

With[{n = size - 1},
   Block[{m1, m2, ones},
    m1 = LowerTriangularize[ConstantArray[1, {n, n}], -1];
    m2 = m1 - DiagonalMatrix[Range[n]];
    ones = Developer`ToPackedArray@{ConstantArray[1, n]};
    ArrayFlatten[{{1, ones}, {Transpose@ones, m2}}]
    ]]; // RepeatedTiming
(*  {0.000041, Null}  *)

Using ArrayPad[m2, {{1, 0}, {1, 0}}, 1], as suggested by @kglr, is about the same speed (the timings fluctuate around each other), even though it eliminates a line of code (for ones). But ones is just 1 x 100 and rather small compared to the matrices. At size = 1000, ArrayFlatten[] took 0.014-0.015 sec., ArrayPad[] took 0.015-0.16 sec., and my method took 0.012-0.014 sec.

My own answer avoids unpacking and uses vectorized functions. Outer[Plus, array1, array2] and Outer[List, array1, array2] are extremely efficient special cases of Outer[] on packed arrays.

Appendix: Timing code

sa[size]; // RepeatedTiming                  (* @kglr *)
Table[Which[i == 1, 1, i > j, 1, i == j, -i + 1, True, 0],
   {i, size}, {j, size}]; // RepeatedTiming  (* @Rohit Namjoshi *)
With[{n = size - 1},
 Block[{m1, m2, ones},
  m1 = LowerTriangularize[ConstantArray[1, {n, n}], -1];
  m2 = m1 - DiagonalMatrix[Range[n]];
  ones = {ConstantArray[1, n]};
  ArrayFlatten[{{1, ones}, {Transpose@ones, m2}}]
  ]]; // RepeatedTiming                      (* @MikeY *)
With[{r    = Range@size,
      zero = Developer`ToPackedArray@{0}},
   UnitStep[Outer[Plus, r, -r] + PadRight[ConstantArray[r, {1}], {size, size}]] - 
    DiagonalMatrix@PadRight[zero, size, r]
   ]; // RepeatedTiming                      (* @Michael E2 *)

You can also use SparseArray as follows

sa[n_] := SparseArray[{{i_, i_} /; i > 1 -> 1 - i, {i_, j_} /; 1 < i < j ->  0}, {n, n}, 1]

matrix2 = sa[7] 
matrix2 == matrix

True

Tags:

Matrix

Recent Posts