Home Page
Site Map
Login/out
Message Board
Download(c) 2005-2025 XLOG Technologies AG
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="UTF-8">
<title editable="nocomment">Example 62</title>
<link rel="stylesheet" href="style.css" type="text/css">
</head>
<!-- Warranty & Liability -->
<!-- To the extent permitted by applicable law and unless explicitly -->
<!-- otherwise agreed upon, XLOG Technologies AG makes no warranties -->
<!-- regarding the provided information. XLOG Technologies AG assumes -->
<!-- no liability that any problems might be solved with the information -->
<!-- provided by XLOG Technologies AG. -->
<!-- -->
<!-- Rights & License -->
<!-- All industrial property rights regarding the information - copyright -->
<!-- and patent rights in particular - are the sole property of XLOG -->
<!-- Technologies AG. If the company was not the originator of some -->
<!-- excerpts, XLOG Technologies AG has at least obtained the right to -->
<!-- reproduce, change and translate the information. -->
<!-- -->
<!-- Reproduction is restricted to the whole unaltered document. -->
<!-- Reproduction of the information is only allowed for non-commercial -->
<!-- uses. Selling, giving away or letting of the execution of the -->
<!-- library is prohibited. The library can be distributed as part of -->
<!-- your applications and libraries for execution provided this comment -->
<!-- remains unchanged. -->
<!-- -->
<!-- Restrictions -->
<!-- Only to be distributed with programs that add significant and primary -->
<!-- functionality to the library. Not to be distributed with additional -->
<!-- software intended to replace any components of the library. -->
<!-- -->
<!-- Trademarks -->
<!-- Jekejeke is a registered trademark of XLOG Technologies AG. -->
<body>
<h1>Example 62: Vector Matrix</h1>
<h2>Introduction</h2>
<p>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.</p>
<h2>Vector Representation</h2>
<p>Vectors will be represented as lists <code>[A1,..,An]</code>:</p>
<ul>
<li>Where <code>A1 .. An</code> are numbers.</li>
</ul>
<p>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:</p>
<pre>
[W1,..,Wn] . [X1,..,Xn] = W1*X1 + .. + Wn*Xn
</pre>
<p>The above index based definition leads to a list based implementation
where we first form the list of products <code>[W1*X1,..,Wn*Xn]</code>
and then sum over the products. The Prolog code reads as follows:</p>
<pre class="code">
:- 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).
</pre>
<p>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 <code>0</code>, and the subsequently form a product and a sum
<code>.. ((0+W1*X1)+W2*X2) .. +Wn*Xn</code>. This leads to an alternative
more direct implementation:</p>
<pre class="code">
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).
</pre>
<h2>Matrice Representation</h2>
<p>Matrixes will be represeted as lists <code>[V1,..,Vm]</code>:</p>
<ul>
<li>Where <code>V1 .. Vm</code> are vectors.</li>
</ul>
<p>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:</p>
<pre>
(A^T)ij = Aji
</pre>
<p>Surprisingly we can use list processing to compute the transpose:</p>
<pre class="code">
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).
</pre>
<p>Linear algebra gets its name from linear equations of the form
<code>W1*X1 + .. + Wn*Xn = B</code>. If they appear in plenty the
weights <code>Wij</code> and targets Bi can be expressed as a matrix
A respectively a vector y. And the equations can be expressed as a
linear mapping <code>A x = y</code>.</p>
<pre class="code">
matapply(M, V, W) :-
maplist(vecdot(V), M, W).
?- matapply([[0,1],[1,0]], [2,3], X).
</pre>
<p>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 <code>A B = C</code>. Again list
processing is enough to implement this operation:</p>
<pre class="code">
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).
</pre>
<h2>Matrix Exponentiation</h2>
<p>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:</p>
<pre>
F0 = 0, F1 = 1, Fn+2 = Fn+1 + Fn
</pre>
<p>We can use a vector <code>Vn = [Fn+1, Fn]</code>. The initial condition
can now be expressed as <code>V0 = [1, 0]</code>. For the recurrence
relation we need to find a matrix so that <code>M * Vn = Vn+1</code>.
The following matrix does that:</p>
<pre>
| 1 1 | | Fn+1 | | Fn+1 + Fn | | Fn+2 |
| | | | = | | = | |
| 1 0 | | Fn | | Fn+1 | | Fn+1 |
</pre>
<p>Lets check whether we can reproduce F10 = 55, F11 = 89 and F12 = 144:</p>
<pre class="code">
?- matapply([[1,1],[1,0]], [89,55], X).
</pre>
<p>In general <code>M^n V0 = Vn</code>, a matrix exponentiation M^n can be implemented:</p>
<pre class="code">
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).
</pre>
<p>We can again reproduce F12 = 144:</p>
<pre class="code">
fib(N, X) :-
matexp([[1,1],[1,0]], N, H),
matapply(H, [1,0], [_,X]).
?- fib(12, X).
</pre>
<p>What are the last 8 digits of F1000000?</p>
<pre class="code">
?- fib(1000000, _X), Y is _X mod 10^8.
</pre>
<p>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.</p>
<h2>Conclusions</h2>
<p>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)?".</p>
<script type="module">
import {
notebook_async
} from "../../../../dogelog/player/canned/dogelog.mjs";
await notebook_async();
</script>
</body>
</html>