10
$\begingroup$

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:

enter image description here

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.

$\endgroup$

5 Answers 5

15
$\begingroup$
size = 7;
Table[Which[i == 1, 1, i > j, 1, i == j, -i + 1, True, 0], {i, size}, {j, size}]
$\endgroup$
2
  • $\begingroup$ Worked very well, now I'll have to understand your syntax =) $\endgroup$
    – Allure
    Jul 2, 2019 at 3:41
  • $\begingroup$ You might be interested in a tweak to your method explained in my answer. $\endgroup$
    – Michael E2
    Jul 2, 2019 at 17:35
12
$\begingroup$

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 *)
$\endgroup$
7
  • $\begingroup$ Nice comparison, thanks. I chose my approach more for ease of understanding, or maybe to align with the underlying logic of the matrix, rather than performance. Nice to see it did alright on the performance front too. I learned something! $\endgroup$
    – MikeY
    Jul 2, 2019 at 17:41
  • $\begingroup$ Thanks Michael. I think you meant "can be improved by using Array". $\endgroup$ Jul 2, 2019 at 17:50
  • $\begingroup$ @RohitNamjoshi Thanks for pointing out the error. $\endgroup$
    – Michael E2
    Jul 2, 2019 at 18:06
  • $\begingroup$ @MikeY You could also use ones = ConstantArray[1, {1, n}] instead of ToPackedArray. $\endgroup$
    – Michael E2
    Jul 2, 2019 at 18:14
  • $\begingroup$ @MichaelE2, I tweaked my answer, thanks for the tip. $\endgroup$
    – MikeY
    Jul 2, 2019 at 19:12
9
$\begingroup$

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

$\endgroup$
9
$\begingroup$

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

$\endgroup$
1
  • 3
    $\begingroup$ you can also use ArrayPad[#, {{1, 0}, {1, 0}}, 1] & on m2 for the last step (+1) $\endgroup$
    – kglr
    Jul 2, 2019 at 14:38
4
$\begingroup$

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}

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.