Our old Prolog system had a little computer algebra system (CAS) that supported vectors and marices. In the following we report on a new take that is less leaning towards object oriented programming and more leaning towards higher order logic programming.
Vectors will be represented as lists [A1,..,An]:
A1 .. An are numbers.If we would use index based access via a list predicate such as nth1/3, this would be slow. We demonstrate how to realize some standard canon of linear algebra operations with higher order list processing predicates only. We first look at the inner product:
[W1,..,Wn] . [X1,..,Xn] = W1*X1 + .. + Wn*Xn
The above index based definition leads to a list based implementation
where we first form the list of products [W1*X1,..,Wn*Xn]
and then sum over the products. The Prolog code reads as follows:
:- ensure_loaded(library(lists)).
vecdot2(V, W, X) :-
maplist(mul, V, W, H),
foldl(add, H, 0, X).
mul(X, Y, Z) :- Z is X*Y.
add(X, Y, Z) :- Z is X+Y.
?- vecdot2([1,3,-5], [4,-2,-1], X).
We might want to avoid the materialization of the list H. It helps
to view the sum in the form of a Horner schema, where we start with
zero 0, and the subsequently form a product and a sum
.. ((0+W1*X1)+W2*X2) .. +Wn*Xn. This leads to an alternative
more direct implementation:
vecdot(V, W, X) :- foldl(muladd, V, W, 0, X). muladd(X, Y, Z, T) :- T is X*Y+Z. ?- vecdot([1,3,-5], [4,-2,-1], X).
Matrixes will be represeted as lists [V1,..,Vm]:
V1 .. Vm are vectors.It is assume that the vectors V1 .. Vm have all the same size n, leading to a m x n matrix. The vectors are called the rows of the matrix, and the vector elements run along the column indexes. An operation that doesn't involve arithmetic is taking the transpose:
(A^T)ij = Aji
Surprisingly we can use list processing to compute the transpose:
mattran([H|T], R) :- mattran(T, H, R). mattran([], U, R) :- maplist(rdot([]), U, R). mattran([J|T], H, R) :- mattran(T, J, TC), maplist(rdot, TC, H, R). rdot(T, H, [H|T]). ?- mattran([[1,2],[3,4],[5,6]], X).
Linear algebra gets its name from linear equations of the form
W1*X1 + .. + Wn*Xn = B. If they appear in plenty the
weights Wij and targets Bi can be expressed as a matrix
A respectively a vector y. And the equations can be expressed as a
linear mapping A x = y.
matapply(M, V, W) :- maplist(vecdot(V), M, W). ?- matapply([[0,1],[1,0]], [2,3], X).
If there is a multiplicity of vectors x and y, we can arrange
them into matrix B and C. The system of equations can then be
expressed as matrix multiplication A B = C. Again list
processing is enough to implement this operation:
matmul(M, N, R) :- mattran(N, H), maplist(matapply(H), M, R). ?- matmul([[1,2,3],[3,2,1]], [[4,5],[6,5],[4,6]], X).
We will turn to an application of linear algbra to number theory. It is possible to expess recurence relation such as the one of Fibonnaci numbers by a matrix. The initial condition and the recurrence relation is:
F0 = 0, F1 = 1, Fn+2 = Fn+1 + Fn
We can use a vector Vn = [Fn+1, Fn]. The initial condition
can now be expressed as V0 = [1, 0]. For the recurrence
relation we need to find a matrix so that M * Vn = Vn+1.
The following matrix does that:
| 1 1 | | Fn+1 | | Fn+1 + Fn | | Fn+2 |
| | | | = | | = | |
| 1 0 | | Fn | | Fn+1 | | Fn+1 |
Lets check whether we can reproduce F10 = 55, F11 = 89 and F12 = 144:
?- matapply([[1,1],[1,0]], [89,55], X).
In general M^n V0 = Vn, a matrix exponentiation M^n can be implemented:
matexp(M, 1, M) :- !. matexp(M, N, R) :- N mod 2 =:= 0, !, I is N // 2, matexp(M, I, H), matmul(H, H, R). matexp(M, N, R) :- I is N-1, matexp(M, I, H), matmul(H, M, R).
We can again reproduce F12 = 144:
fib(N, X) :-
matexp([[1,1],[1,0]], N, H),
matapply(H, [1,0], [_,X]).
?- fib(12, X).
What are the last 8 digits of F1000000?
?- fib(1000000, _X), Y is _X mod 10^8.
The above will have a huge intermediate result _X, which the top-level doesn't show since we used a marked variable. This means we used bigint computation and only took modulo at the end. Alternatively we could have used modulo computation during the matrix exponentiation already.
Currently we only support numerical computation, either floats or bigint. We could implement Prolog operations that have strict vector or matrix signatures. Running it in the browser, we can answer in a blink "What are the last 8 digits of fib(1000000)?".