Matrix Arithmetic

iThink and STELLA support a variety of matrix operations. In this post, I will explain the following common operations and how to perform them:

• Addition and Subtraction
• Transposition
• Extracting the Diagonal of a Square Matrix
• Multiplying and Dividing by a Scalar
• Multiplying a Square Matrix by a Vector
• Multiplying a Non-Square Matrix by a Vector
• Multiplying Two Matrices

In the examples provided, I use five different arrays of varying dimensions. They are all defined using the subscript name Dim3, which is size 3.

a: 3×3
b: 3×3
c: 3×3
d: 3×1 [one-dimensional]
e: 3×1 [one-dimensional]

I also use two arrays based on Dim3, as well as the subscript name Dim2, which is size 2.

s: 2×3
t: 3×2

Addition and Subtraction

The easiest functions to perform are addition and subtraction. When two arrays have the same dimensions, they can be added and subtracted in the usual way, placing the result in a third array of the same dimensions.

For example, to place the difference of b and c into array a, set the equation of a to

b[Dim3, Dim3] – c[Dim3, Dim3]

Transposition

Transposition of a 2-dimensional array is also relatively easy, especially when the dimensions differ. Suppose we wish to transpose s and place the result in t. Simply set the equation for t to

s[Dim2, Dim3]

This will properly fill t, which has dimensions Dim3xDim2, with:

s₁₁   s₂₁
s₁₂   s₂₂
s₁₃   s₂₃

However, if the two dimensions of the array have the same name, the ARRAYVALUE built-in must be used. Suppose we wish to transpose b, placing its transposition in a. Set the equation of a to

ARRAYVALUE(b[*, *], ARRAYIDX(2), ARRAYIDX())

Note the index of a’s second dimension is used as the index for b’s first dimension and the index of a’s first dimension is used as the index for b’s second dimension.

Extracting the Diagonal of a Square Matrix

Extracting the diagonal from a 2-dimensional array is also very easy. To place the diagonal of a into d, set the equation of d to

a[Dim3, Dim3]

Multiplying and Dividing by a Scalar

Multiplying or dividing by a scalar works as expected. The contents of b can be multiplied by a scalar x, for example, by setting the equation of a to

x*b[Dim3, Dim3]

Multiplying a Square Matrix by a Vector

Matrix multiplication, as defined in linear algebra, is not directly supported at present. It is possible to multiply two arrays, but the results will not match matrix multiplication. When you multiply a 1-dimensional array by a 2-dimensional array, the size of the 1-dimensional array must match either the number of rows or the number of columns in the 2-dimensional array. In both cases, iThink and STELLA treats the 1-dimensional array as either a column or row of scalars that are multiplied across each corresponding row or column, respectively, of the 2-dimensional array.

Given the following equation in a,

d[Dim3]*b[Dim3, Dim3]

d will be treated as a column of scalars (when the 2-dimensional array is square with the same dimension name in both directions, the 1-dimensional array is always treated as a column vector). Each scalar in d will be multiplied across each corresponding row of b. Note this equation does not represent a valid matrix multiplication, and also that switching the order to b[Dim3, Dim3]*d[Dim3], while now a valid matrix multiplication, will not change the result. I recommend leaving the order as above to emphasize that it is not a normal matrix multiplication.

More precisely, the array a will be populated with these values:

d₁*b₁₁   d₁*b₁₂   d₁*b₁₃
d₂*b₂₁   d₂*b₂₂   d₂*b₂₃
d₃*b₃₁   d₃*b₃₂   d₃*b₃₃

Note that these are not the correct terms for matrix multiplication. If we multiply b[Dim3, Dim3]*d[Dim3], which is the only correct way to express it, the linear algebra result is a one-dimensional vector with these values:

d₁*b₁₁ + d₂*b₁₂ + d₃*b₁₃
d₁*b₂₁ + d₂*b₂₂ + d₃*b₂₃
d₁*b₃₁ + d₂*b₃₂ + d₃*b₃₃

If we desire to generate these terms, we must first transpose b by placing the transposition equation given above in c:

ARRAYVALUE(b[*, *], ARRAYIDX(2), ARRAYIDX())

This puts the following values in c, i.e., b transposed:

b₁₁   b₂₁   b₃₁
b₁₂   b₂₂   b₃₂
b₁₃   b₂₃   b₃₃

The equation for a can now be set to

d[Dim3]*c[Dim3, Dim3]

This equation will populate a with these values:

d₁*b₁₁   d₁*b₂₁   d₁*b₃₁
d₂*b₁₂   d₂*b₂₂   d₂*b₃₂
d₃*b₁₃   d₃*b₂₃    d₃*b₃₃

Note that this gives the proper terms for the matrix multiplication, arranged vertically. The final answer is found by summing each column, using the following equation in vector e:

ARRAYSUM(a[*, Dim3])

Aside: The terms we get from the original equation for a, d[Dim3]*b[Dim3, Dim3], are the correct terms for d^Tb, i.e., d-transposed times b, again arranged vertically. If this was the desired result, we do not need to first transpose b and can just sum the columns into vector e as shown above.

Multiplying a Non-Square Matrix by a Vector

If the matrices do not use the same dimension name for both dimensions (square or not), you have much more control over the resultant terms and do not have to resort to transposition at all. The following discussion uses the following additional arrays:

f: 2×3
g: 2×1 [one-dimensional]
s: 2×3t: 3×2

If we multiply f times d, we will get a different result based on whether the resultant array is 2×3 or 3×2. Placing the equation d[Dim3]*f [Dim2, Dim3] into array s yields exactly the terms needed for a matrix multiplication, arranged horizontally:

d₁*f₁₁   d₂*f₁₂   d₃*f₁₃
d₁*f₂₁   d₂*f₂₂   d₃*f₂₃

Setting the equation of g to

ARRAYSUM(s[Dim2, *])

will properly sum the rows into g, giving the result of the matrix multiplication f *d.

If, instead, we place the same equation d[Dim3]*f [Dim2, Dim3] into array t, we get these terms instead:

d₁*f₁₁   d₁*f₂₁
d₂*f₁₂   d₂*f₂₂
d₃*f₁₃   d₃*f₂₃

The correct terms are still available, only now they are arranged vertically and the columns must be summed to arrive at the proper result:

ARRAYSUM(t[*, Dim2])

Multiplying by g instead, i.e., placing the equation g[Dim2]*f [Dim2, Dim3] into s, yields these terms:

g₁*f₁₁   g₁*f₁₂   g₁*f₁₃
g₂*f₂₁   g₂*f₂₂   g₂*f₂₃

These are the terms for g^Tf , arranged vertically. Placing the equation into t gives the same terms horizontally. To calculate f ^Tg, it is necessary to first transpose f by setting t to the equation

f [Dim2, Dim3]

When this is multiplied by g, the correct terms for f ^Tg will be generated (horizontally in a 3×2 array and vertically in a 2×3 array).

Multiplying Two Matrices

Unfortunately, there is no way to easily generate the terms for multiplying two matrices. The equations must be individually tailored for each case. As an example, let’s look at multiplying t*s and storing the result in a. Array t is 3×2 and array s is 2×3, so the result will fit in a (3×3). The correct answer is:

t₁₁*s₁₁ + t₁₂*s₂₁    t₁₁*s₁₂ + t₁₂*s₂₂    t₁₁*s₁₃ + t₁₂*s₂₃
t₂₁*s₁₁ + t₂₂*s₂₁    t₂₁*s₁₂ + t₂₂*s₂₂    t₂₁*s₁₃ + t₂₂*s₂₃
t₃₁*s₁₁ + t₃₂*s₂₁    t₃₁*s₁₂ + t₃₂*s₂₂    t₃₁*s₁₃ + t₃₂*s₂₃

Examining these values for a pattern that could be used in an apply-to-all equation, we can see the proper equation for array a is

t[Dim3, 1]*s[1, Dim3] + t[Dim3, 2]*s[2, Dim3]

This pattern, incidentally, will hold in general for any two matrices. However, this equation does not work as an apply-to-all equation because the same name is used for both the first and second dimensions so it is not possible to tell which of the two varying Dim3 values to use. In this case, iThink and STELLA currently always substitutes the value of the first dimension, which will not give the correct answer. This will always be a problem when the resulting array is square. One workaround is to define an additional dimension name that is the same size as Dim3. This works, but may lead to other array incompatibilities. The other solution is similar to how we transposed a square matrix – explicitly use ARRAYVALUE for any references to the second dimension:

t[Dim3, 1]*ARRAYVALUE(s[*, *], 1, ARRAYIDX(2))
+ t[Dim3, 2]*ARRAYVALUE(s[*, *], 2, ARRAYIDX(2))

This is certainly not as clear, but it works and is only necessary when the result is square. One last example underscores this point and also shows how to expand the general pattern to a larger array. For this example, one more dimension name, Dim4, of size 4 and two more arrays are needed:

u: 3×4
v: 2×4

To calculate the result of s*u in array v, where s is 2×3, u is 3×4, and v is (correctly) 2×4, set the equation for v to:

s[Dim2, 1]*u[1, Dim4] + s[Dim2, 2]*u[2, Dim4] + s[Dim2, 3]*u[3, Dim4]

Note this same method can be used for multiplying a 2-dimensional array by a 1-dimensional array instead of the method presented above. It avoids both the intermediate transposition and product arrays, but results in more complex expressions that must be tailored to each combination of array sizes.