```    1:- module(linear_algebra,
2          [   matrix_dimensions/3,      % ?Matrix, ?Rows, ?Columns
3              matrix_identity/2,        % +Order, -Matrix
4              matrix_transpose/2,       % ?Matrix0, ?Matrix
5              matrix_multiply_vector/3,
6              matrix_multiply/3,
7              matrix_rotation/2,        % ?Theta, ?Matrix
8
9              vector_dimension/2,
10              vector_distance/2,        % ?V, ?Distance
11              vector_distance/3,        % ?U, ?V, ?Distance
12              vector_translate/3,       % ?U, ?V, ?W
13              vector_multiply/3,
14              vector_scale/3,           % ?Scalar, ?U, ?V
16
17              scalar_translate/3,
18              scalar_multiply/3,
19              scalar_power/3            % ?X, ?Y, ?Z
20          ]).   21
22:- use_module(library(clpr)).   23:- use_module(library(clpq), []).```

# Linear algebra

"The introduction of numbers as coordinates is an act of violence."--Hermann Weyl, 1885-1955.

Vectors are just lists of numbers, or scalars. These scalars apply to arbitrary abstract dimensions. For example, a two-dimensional vector [1, 2] applies two scalars, 1 and 2, to dimensional units i and j; known as the basis vectors for the coordinate system.

Is it possible, advisable, sensible to describe vector and matrix operations using Constraint Logic Programming (CLP) techniques? That is, since vectors and matrices are basically columns and rows of real-numeric scalars, their operators amount to constrained relationships between real numbers and hence open to the application of CLP over reals. The simple answer is yes, the linear_algebra predicates let you express vector operators using real-number constraints.

Constraint logic adds some important features to vector operations. Suppose for instance that you have a simple addition of two vectors, a vector translation of U+V=W. Add U to V giving W. The following statements all hold true. Note that the CLP-based translation unifies correctly when W is unknown but also when U or V is unknown. Given any two, you can ask for the missing vector.

```?- vector_translate([1, 1], [2, 2], W).
W = [3.0, 3.0] ;
false.
?- vector_translate([1, 1], V, [3, 3]).
V = [2.0, 2.0] ;
false.
?- vector_translate(U, [2, 2], [3, 3]).
U = [1.0, 1.0] ;
false.```

Note also that the predicate answers non-deterministically with back-tracking until no alternative answer exists. This presumes that alternatives could exist at least in theory if not in practice. Trailing choice-points remain unless you cut them.

/

matrix_dimensions(?Matrix:list(list(number)), ?Rows:nonneg, ?Columns:nonneg) is semidet
Dimensions of Matrix where dimensions are Rows and Columns.

A matrix of M rows and N columns is an M-by-N matrix. A matrix with a single row is a row vector; one with a single column is a column vector. Because the linear_algebra module uses lists to represent vectors and matrices, you need never distinguish between row and column vectors.

Boundary cases exist. The dimensions of an empty matrix [] equals [0, _] rather than [0, 0]. And this works in reverse; the matrix unifying with dimensions [0, _] equals [].

```   83matrix_dimensions(Matrix, Rows, Columns) :-
84    length(Matrix, Rows),
85    matrix_dimensions_(Matrix, Columns).
86
87matrix_dimensions_([], _).
88matrix_dimensions_([V0|V], Columns) :-
89    vector_dimension(V0, Columns),
90    matrix_dimensions_(V, Columns).```
matrix_identity(+Order:nonneg, -Matrix:list(list(number))) is semidet
Matrix becomes an identity matrix of Order dimensions. The result is a square diagonal matrix of Order rows and Order columns.

The first list of scalars (call it a row or column) becomes 1 followed by Order-1 zeros. Subsequent scalar elements become an Order-1 identity matrix with a 0-scalar prefix for every sub-list. Operates recursively albeit without tail recursion.

Fails when matrix size Order is less than zero.

```  105matrix_identity(0, []) :-
106    !.
107matrix_identity(Order, [[1|Vector]|Matrix]) :-
108    integer(Order),
109    Order > 0,
110    Order0 is Order - 1,
111    vector_dimension(Vector, Order0),
112    maplist(=(0), Vector),
113    matrix_identity(Order0, Matrix0),
114    maplist([Scalars, [0|Scalars]]>>true, Matrix0, Matrix).```
matrix_transpose(?Matrix0:list(list(number)), ?Matrix:list(list(number))) is semidet
Transposes matrices. The matrix is a list of lists. Fails unless all the sub-lists share the same length. Works in both directions, and works with non-numerical elements. Only operates at the level of two-dimensional lists, a list with sub-lists. Sub-sub-lists remain lists and un-transposed if sub-lists comprise list elements.
```  125matrix_transpose([U0|U], []) :-
126    maplist(=([]), [U0|U]),
127    !.
128matrix_transpose([U0|U], [V0|V]) :-
129    maplist([[H|T], H, T]>>true, [U0|U], V0, U_),
130    matrix_transpose(U_, V).
131
132matrix_multiply_vector(M, U, V) :-
133    maplist(vector_scale, U, M, [W0|W]),
134    foldl(vector_translate, W, W0, V).
135
136matrix_multiply(A, B, M) :-
137    maplist(matrix_multiply_vector(B), A, M).```
matrix_rotation(?Theta:number, ?Matrix:list(list(number))) is nondet
The constructed matrix applies to column vectors [X, Y] where positive Theta rotates X and Y anticlockwise; negative rotates clockwise. Transpose the rotation matrix to reverse the angle of rotation; positive for clockwise, negative anticlockwise.
```  147matrix_rotation(Theta, [[A, B], [C, A]]) :-
148    {A =:= cos(Theta), B =:= sin(Theta), C =:= -B}.
149
150vector_dimension(V, Columns) :- length(V, Columns).```
vector_distance(?V:list(number), ?Distance:number) is semidet
vector_distance(?U:list(number), ?V:list(number), ?Distance:number) is semidet
Distance of the vector V from its origin. Distance is Euclidean distance between two vectors where the first vector is the origin. Note that Euclidean is just one of many distances, including Manhattan and chessboard, etc. The predicate is called distance, rather than length. The term length overloads on the dimension of a vector, its number of numeric elements.
```  163vector_distance(V, Distance) :-
164    maplist(scalar_power(2), V, [X0|X]),
165    foldl(scalar_translate, X, X0, Scalar),
166    scalar_power(2, Distance, Scalar),
167    !.
168
169vector_distance(U, V, Distance) :-
170    vector_translate(U, W, V),
171    vector_distance(W, Distance).```
vector_translate(?U, ?V, ?W) is nondet
Translation works forwards and backwards. Since U+V=W it follows that U=W-V and also V=W-U. So for unbound U, the vector becomes W-V and similarly for V.
```  179vector_translate([], [], []).
180vector_translate([X|U], [Y|V], [Z|W]) :-
181    scalar_translate(X, Y, Z),
182    vector_translate(U, V, W).
183
184vector_multiply([], [], []).
185vector_multiply([X|U], [Y|V], [Z|W]) :-
186    scalar_multiply(X, Y, Z),
187    vector_multiply(U, V, W).```
vector_scale(?Scalar:number, ?U:list(number), ?V:list(number)) is nondet
Vector U scales by Scalar to V.

What is the difference between multiply and scale? Multiplication multiplies two vectors whereas scaling multiplies a vector by a scalar; hence the verb to scale. Why is the scalar at the front of the argument list? This allows the meta-call of `vector_scale(Scalar)` passing two vector arguments, e.g. when mapping lists of vectors.

The implementation performs non-deterministically because the CLP(R) library leaves a choice point when searching for alternative arithmetical solutions.

```  204vector_scale(_, [], []).
205vector_scale(Scalar, [X|U], [Y|V]) :-
206    scalar_multiply(Scalar, X, Y),
207    vector_scale(Scalar, U, V).```
Heading in radians of vector V. Succeeds only for two-dimensional vectors. Normalises the Heading angle in (+, -) mode; negative angles wrap to the range between pi and two-pi. Similarly, normalises the vector V in (-, +) mode; V has unit length.
```  216vector_heading([X, Y], Heading) :-
218    !,
219    Angle is atan2(Y, X),
220    (   Angle < 0
221    ->  Heading is Angle + 2 * pi
223    ).
The first argument X is the exponent rather than Y, first rather than second argument. This allows you to curry the predicate by fixing the first exponent argument. In other words, ```scalar_power(2, A, B)``` squares A to B.
`  241scalar_power(X, Y, Z) :- {Z =:= pow(Y, X)}`