## 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

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 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 *Dim3*x*Dim2*, with:

*s*_{11} *s*_{21}

*s*_{12} *s*_{22}

*s*_{13} *s*_{23}

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*_{1}**b*_{11} *d*_{1}**b*_{12} *d*_{1}**b*_{13}

*d*_{2}**b*_{21} *d*_{2}**b*_{22} *d*_{2}**b*_{23}

*d*_{3}**b*_{31} *d*_{3}**b*_{32} *d*_{3}**b*_{33}

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*_{1}**b*_{11} + *d*_{2}**b*_{12} + *d*_{3}**b*_{13}

*d*_{1}**b*_{21} + *d*_{2}**b*_{22} + *d*_{3}**b*_{23}

*d*_{1}**b*_{31} + *d*_{2}**b*_{32} + *d*_{3}**b*_{33}

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*_{11} *b*_{21} *b*_{31}

*b*_{12} *b*_{22} *b*_{32}

*b*_{13} *b*_{23} *b*_{33}

The equation for *a* can now be set to

*d*[*Dim3*]**c*[*Dim3*, *Dim3*]

This equation will populate *a* with these values:

*d*_{1}**b*_{11} *d*_{1}**b*_{21} *d*_{1}**b*_{31}

*d*_{2}**b*_{12} *d*_{2}**b*_{22} *d*_{2}**b*_{32}

*d*_{3}**b*_{13} *d*_{3}**b*_{23} *d*_{3}**b*_{33}

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*^{T}*b*, 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×3*t*: 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*_{1}**f*_{11} *d*_{2}**f*_{12} *d*_{3}**f*_{13}

*d*_{1}**f*_{21} *d*_{2}**f*_{22} *d*_{3}**f*_{23}

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*_{1}**f*_{11} *d*_{1}**f*_{21}

*d*_{2}**f*_{12} *d*_{2}**f*_{22}

*d*_{3}**f*_{13} *d*_{3}**f*_{23}

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*_{1}**f*_{11} *g*_{1}**f*_{12} *g*_{1}**f*_{13}

*g*_{2}**f*_{21} *g*_{2}**f*_{22} *g*_{2}**f*_{23}

These are the terms for *g*^{T}*f* , arranged vertically. Placing the equation into *t* gives the same terms horizontally. To calculate *f* ^{T}*g*, 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* ^{T}*g* will be generated (horizontally in a 3×2 array and vertically in a 2×3 array).** **

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*_{11}**s*_{11} + *t*_{12}**s*_{21} *t*_{11}**s*_{12} + *t*_{12}**s*_{22} *t*_{11}**s*_{13} + *t*_{12}**s*_{23}

*t*_{21}**s*_{11} + *t*_{22}**s*_{21} *t*_{21}**s*_{12} + *t*_{22}**s*_{22} *t*_{21}**s*_{13} + *t*_{22}**s*_{23}

*t*_{31}**s*_{11} + *t*_{32}**s*_{21} *t*_{31}**s*_{12} + *t*_{32}**s*_{22} *t*_{31}**s*_{13} + *t*_{32}**s*_{23}

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.

*If you enjoyed this post, make sure you subscribe to my RSS feed!*