# Course:CPSC312-2018-Matrixlog

## Matrixlog: Programming with Matrices and Logic

#### Authors:

Jerry Yin, James Johnson

### What is the problem?

A matrix is a multidimensional array of numbers. Matrices are commonly used in scientific computations. One application of matrices is that they are used to represent linear systems of equations. For example the linear system of equations with ${\displaystyle n}$ unknowns (${\displaystyle x_{i}}$)

${\displaystyle {\begin{array}{lcl}a_{1,1}x_{1}+a_{1,2}x_{2}+\ldots +a_{1,n}x_{n}&=&b_{1}\\a_{2,1}x_{1}+a_{2,2}x_{2}+\ldots +a_{2,n}x_{n}&=&b_{2}\\&\vdots \\a_{n,1}x_{1}+a_{n,2}x_{2}+\ldots +a_{n,n}x_{n}&=&b_{n}\end{array}}}$
can be represented as the matrix equation

${\displaystyle {\begin{bmatrix}a_{1,1}&a_{1,2}&\cdots &a_{1,n}\\a_{2,1}&a_{2,2}&\cdots &a_{2,n}\\\vdots &\vdots &\ddots &\vdots \\a_{n,1}&a_{n,2}&\cdots &a_{n,n}\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{2}\\\vdots \\x_{n}\end{bmatrix}}={\begin{bmatrix}b_{1}\\b_{2}\\\vdots \\b_{n}\end{bmatrix}}.}$

We will implement various 2D matrix operations in Prolog. This will include:

• matrix-scalar multiplication,
• matrix-vector multiplication,
• matrix-matrix multiplication,
• matrix transpose, and
• at least one matrix solver.

For the matrix solver, we will implement LU decomposition, a method that decomposes a square invertible matrix ${\displaystyle A}$ into an upper triangular matrix ${\displaystyle U}$ and a lower triangular matrix ${\displaystyle L}$.

${\displaystyle A=LU={\begin{bmatrix}1&0&0&\cdots &0\\l_{2,1}&1&0&\cdots &0\\l_{3,1}&l_{3,2}&1&\cdots &0\\\vdots &\vdots &\vdots &\ddots &\vdots \\l_{n,1}&l_{n,2}&l_{n,3}&\cdots &1\end{bmatrix}}{\begin{bmatrix}u_{1,1}&u_{1,2}&u_{1,3}&\cdots &u_{1,n}\\0&u_{2,2}&u_{2,3}&\cdots &u_{2,n}\\0&0&u_{3,3}&\cdots &u_{3,n}\\\vdots &\vdots &\vdots &\ddots &\vdots \\0&0&0&\cdots &u_{n,n}\end{bmatrix}}}$

After performing this decomposition, we can observe that

{\displaystyle {\begin{aligned}A\mathbf {x} &=\mathbf {b} \\LU\mathbf {x} &=\mathbf {b} \\U\mathbf {x} &=L^{-1}\mathbf {b} \\\mathbf {x} &=U^{-1}(L^{-1}\mathbf {b} ).\\\end{aligned}}}

${\displaystyle L^{-1}\mathbf {b} }$ can be obtained using forward substitution (no matrix inversion necessary), and ${\displaystyle U^{-1}(L^{-1}\mathbf {b} )}$ can be obtained using backwards substitution, thus solving the original system ${\displaystyle A\mathbf {x} =\mathbf {b} }$.

For simplicity, we only implement LU decomposition without partial pivoting. Partial pivoting makes the decomposition more numerically stable by exchanging rows, and allows us to solve systems with zeros along the diagonal, such as

${\displaystyle {\begin{bmatrix}0&1\\1&0\end{bmatrix}}\mathbf {x} ={\begin{bmatrix}3\\4\end{bmatrix}}.}$

Since we are only doing a feasibility study, we won't worry about pathological cases such as these.

### What is the something extra?

We investigated using arg/3 and setarg/3 and CLP(FD) for interesting purposes (see below).

### What did we learn from doing this?

##### arg / setarg

Most matrix algorithms require the ability to access the elements of matrices in arbitrary order. This makes Prolog lists (which are effectively linked lists) less suitable for our purposes. Instead we represent the elements of a matrix as a Prolog term. We represent a matrix as

matrix(M = NumRows, N = NumCols,
elements(E11, E12, ..., E1N,
E21, E22, ..., E2N,
⋮     ⋮    ⋱    ⋮
EM1, EM2, ..., EMN))


We can get constant time random access to the elements in the elements term by using arg/3. arg(Index, Term, Value) is true if and only if the Index-th element of Term has value Value.

?- arg(2, t(10,20,30), V).
V = 20.


Unlike Haskell, Prolog supports mutation. Thus Prolog also provides setarg/3, which forces the Index-th element of Term to have value Value (if Index is less than or equal to the arity of Term and greater than or equal to 1), even if it had a different value previously.

?- X = t(1, 2, 3), setarg(1, X, 7).
X = t(7, 2, 3).


Since some algorithms are more convenient over lists, there is a way to convert between them, with the (=..)/2 operator:

?- elements(1, 2, 3, 4, 5) =.. X.
X = [elements, 1, 2, 3, 4, 5].

##### CLP(FD)

CLP(FD) is a library for performing constraint logic programming over finite domains. CLP(FD) provides the (#=)/2 operator, which is similar to is, except that it works in both directions. (In other words, both sides can be expressions.) This also means that Prolog can determine the solution to implicit equations:

?- 12 #= Y + 7.
Y = 5.


Unfortunately, it only works for integers, whereas we support floating point numbers. So we only use (#=) for index computations.

Another shortcoming of CLP(FD) is that it can fail to solve even very simple systems of equations (which is part of why we wrote Matrixlog). For example, it gives up on the following system:

?- A + B #= 2, A + 2 * B #= 3.
A+2*B#=3,
A+B#=2.


and needs a little hint to find the answer:

?- A + B #= 2, A + 2 * B #= 3, A in -100..100.
A = B, B = 1.


However, to solve a system such as

?- 607*A + 492*B + 283*C + 173*D + 928*E + 262*F + 525*G + 548*H + 375*I +  61*J #= 21388,
|  145*A + 475*B + 109*C +  87*D + 348*E + 177*F + 652*G + 584*H + 970*I + 502*J #= 27558,
|  589*A + 118*B + 852*C + 194*D + 773*E + 361*F + 965*G + 520*H + 183*I + 130*J #= 24050,
|  251*A + 432*B + 250*C + 650*D + 711*E + 988*F + 814*G + 840*H + 228*I + 620*J #= 34618,
|  347*A + 329*B + 666*C + 203*D + 106*E + 800*F + 335*G + 686*H + 967*I + 483*J #= 30511,
|  208*A + 378*B + 912*C + 879*D + 613*E + 667*F + 983*G + 192*H +  11*I + 672*J #= 29519,
|  656*A + 610*B + 651*C + 316*D +  16*E + 437*F + 299*G + 613*H + 983*I + 552*J #= 29159,
|  410*A + 874*B +  78*C + 313*D + 734*E + 213*F + 698*G + 927*H + 217*I + 576*J #= 28607,
|  816*A + 928*B + 200*C + 203*D + 213*E + 650*F + 435*G + 851*H + 491*I + 222*J #= 25541,
|  125*A + 381*B + 305*C + 326*D +  49*E + 959*F + 931*G + 392*H + 198*I + 163*J #= 22170,
|  A in 1..10,
|  B in 1..10,
|  C in 1..10,
|  D in 1..10,
|  E in 1..10,
|  F in 1..10,
|  G in 1..10,
|  H in 1..10,
|  I in 1..10,
|  J in 1..10.
A in 1..10,
... + ... + ... * ... + 326*D+49*E+959*F+931*G+392*H+198*I+163*J#=22170,
... + ... + ... * ... + 203*D+213*E+650*F+435*G+851*H+491*I+222*J#=25541,
... + ... + ... * ... + 313*D+734*E+213*F+698*G+927*H+217*I+576*J#=28607,
... + ... + ... * ... + 316*D+16*E+437*F+299*G+613*H+983*I+552*J#=29159,
... + ... + ... * ... + 879*D+613*E+667*F+983*G+192*H+11*I+672*J#=29519,
... + ... + ... * ... + 203*D+106*E+800*F+335*G+686*H+967*I+483*J#=30511,
... + ... + ... * ... + 650*D+711*E+988*F+814*G+840*H+228*I+620*J#=34618,
... + ... + ... * ... + 194*D+773*E+361*F+965*G+520*H+183*I+130*J#=24050,
... + ... + ... * ... + 87*D+348*E+177*F+652*G+584*H+970*I+502*J#=27558,
... + ... + ... * ... + 173*D+928*E+262*F+525*G+548*H+375*I+61*J#=21388,
B in 1..10,
C in 1..10,
D in 1..10,
E in 1..10,
F in 1..10,
G in 1..10,
H in 1..10,
I in 1..10,
J in 1..10.


even giving hints about the ranges doesn't help. (The solution is A = 1, B = 2, C = 3, …)

Our module can solve such systems (although we use inexact machine arithmetic).

?- solve(matrix(10, 10, elements(607, 492, 283, 173, 928, 262, 525, 548, 375,  61,
|                                145, 475, 109,  87, 348, 177, 652, 584, 970, 502,
|                                589, 118, 852, 194, 773, 361, 965, 520, 183, 130,
|                                251, 432, 250, 650, 711, 988, 814, 840, 228, 620,
|                                347, 329, 666, 203, 106, 800, 335, 686, 967, 483,
|                                208, 378, 912, 879, 613, 667, 983, 192,  11, 672,
|                                656, 610, 651, 316,  16, 437, 299, 613, 983, 552,
|                                410, 874,  78, 313, 734, 213, 698, 927, 217, 576,
|                                816, 928, 200, 203, 213, 650, 435, 851, 491, 222,
|                                125, 381, 305, 326,  49, 959, 931, 392, 198, 163)),
|        X,
|        matrix(10, 1, elements(21388, 27558, 24050, 34618, 30511, 29519, 29159,
|                               28607, 25541, 22170))).
X = matrix(10, 1, elements(1.000000000000012, 2.0000000000000115, 2.999999999999996, 3.9999999999999662, 5.0000000000000115, 6.000000000000019, 7.000000000000004, 7.999999999999962, 8.999999999999995, 10.00000000000003)).

We encountered some annoyances when writing numerical algorithms in Prolog. For example, the fact that the is operator is required to evaluate mathematical expressions leads to verbose code, especially if is is called specifically to evaluate an expression whose result is only used once as the input to some other predicate.