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 $n$ unknowns ( $x_{i}$ )

{\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

{\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 addition,
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 $A$ into an upper triangular matrix $U$ and a lower triangular matrix $L$ .

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

{\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}}

$L^{-1}\mathbf {b}$ can be obtained using forward substitution (no matrix inversion necessary), and $U^{-1}(L^{-1}\mathbf {b} )$ can be obtained using backwards substitution, thus solving the original system $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

{\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)).

Disadvantages of Prolog

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.

The lack of a (#=) operator for rationals / floating point numbers also means that our predicates are less general and powerful than they could be. For example, our LU decomposition predicate differentiates between input and output parameters—you supply A and get back L and U—instead of being able to go in "both directions" (supply L and U to get A) or automatically get all solutions (supply just L to get possible U's and A's).

The largest disadvantage, however, is that the syntax of Prolog is more optimized towards declarative queries. This works well when your computation can be efficiently expressed as a declarative query, but as shown in the CLP(FD) section, simply providing enough facts to solve the problem is often not enough to lead Prolog to the answer. In the case of solving systems of equations, it is probably unreasonable to expect Prolog to figure out how to do LU decompositions and so on on its own, so we need to specify the steps of the LU decomposition algorithm. The steps of these matrix algebra algorithms are traditionally presented imperatively (especially within a pedagogical context, so arguably it is more natural to express such algorithms imperatively since that is how most people writing matrix solvers have been taught to think about matrix algorithms), thus trying to re-express these algorithms in a declarative style can result in a lot of friction for the end-user programmer.

Conclusion

In conclusion, we have shown that it is possible to implement traditional matrix algorithms such as matrix multiplication and LU decomposition in Prolog, although our efforts to do so were painful. When implementing traditional matrix algorithms in Prolog, they must be translated into Prolog's declarative style. Implementing these algorithms is necessary, since even when Prolog is given enough information to solve a system of equations, it is often not able to do so without substantial additional guidance from the programmer—it cannot figure out the algorithm to solve systems on its own.

Links to code

https://github.com/enjmiah/matrixlog