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 15 vector_heading/2, % ?V, ?Heading 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), []).
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).
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).
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).
147matrix_rotation(Theta, [[A, B], [C, A]]) :- 148 {A =:= cos(Theta), B =:= sin(Theta), C =:= -B}. 149 150vector_dimension(V, Columns) :- length(V, Columns).
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).
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).
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).
216vector_heading([X, Y], Heading) :- 217 var(Heading), 218 !, 219 Angle is atan2(Y, X), 220 ( Angle < 0 221 -> Heading is Angle + 2 * pi 222 ; Heading = Angle 223 ). 224vector_heading([X, Y], Heading) :- 225 Y is sin(Heading), 226 X is cos(Heading). 227 228scalar_translate(X, Y, Z) :- {Z =:= X + Y}. 229 230scalar_multiply(X, Y, Z) :- {Z =:= X * Y}.
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)}
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.
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.
/