How can I define a very large matrix efficiently?

Question

Here's a 7x7 matrix:

(matrix = {{1, 1, 1, 1, 1, 1, 1}, {1, -1, 0, 0, 0, 0, 0}, {1, 
     1, -2, 0, 0, 0, 0}, {1, 1, 1, -3, 0, 0, 0}, {1, 1, 1, 1, -4, 0, 
     0}, {1, 1, 1, 1, 1, -5, 0}, {1, 1, 1, 1, 1, 
     1, -6}}) // MatrixForm 

The matrix looks like this:

How can I define the 100x100 equivalent of this matrix efficiently without having to manually type out every element? The documentation gives some ideas, but although there's certainly an order in this matrix, it's not something that translates easily (at all?) to the commands used in the documentation.

Answer 1

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

Answer 2

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

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 *)

Answer 3

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

Answer 4

One more, for fun and education. Treating it like a block matrix, do the lower right first. This makes a lower triangular matrix of 1s, with the -1 in the command telling it to make the diagonal = zero.

m1=LowerTriangularize[ConstantArray[1, {6, 6}], -1]

$$ \left( \begin{array}{cccccc} 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 & 0 & 0 \\ 1 & 1 & 1 & 1 & 1 & 0 \\ \end{array} \right) $$

m2 = m1 - DiagonalMatrix[Range[6]]

$$ \left( \begin{array}{cccccc} -1 & 0 & 0 & 0 & 0 & 0 \\ 1 & -2 & 0 & 0 & 0 & 0 \\ 1 & 1 & -3 & 0 & 0 & 0 \\ 1 & 1 & 1 & -4 & 0 & 0 \\ 1 & 1 & 1 & 1 & -5 & 0 \\ 1 & 1 & 1 & 1 & 1 & -6 \\ \end{array} \right) $$

Now the upper right, a $1 \times 6$ matrix of ones

ones = ConstantArray[1, {1, n}]

$$ \left( \begin{array}{cccccc} 1 & 1 & 1 & 1 & 1 & 1 \\ \end{array} \right) $$

Use ArrayFlatten to get the final answer.

ArrayFlatten[{{1, ones}, {Transpose@ones, m2}}]

$$ \left( \begin{array}{ccccccc} 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & -1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & -2 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & -3 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 & -4 & 0 & 0 \\ 1 & 1 & 1 & 1 & 1 & -5 & 0 \\ 1 & 1 & 1 & 1 & 1 & 1 & -6 \\ \end{array} \right) $$

Answer 5

As already suggested in the other answers, vectorization is a key to high performance when working with large numeric arrays. To further increase calculation speed it is useful to remove unnecessary intermediate steps, because memory operations on large arrays can be quite time consuming. Here is a vectorized function with very few intermediate steps:

mFast[n_] := LowerTriangularize[DiagonalMatrix[-Range[n]] + 1] + UnitVector[n, 1];

Timing measurement:

RepeatedTiming[mFast[100];]

{0.000032, Null}

This can be compared to MichaelE2's function:

mMichaelE2[size_] := 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[mMichaelE2[100];]

{0.000047, Null}