1:- module(matrix,
2 [matrix_multiply/3,
3 matrix_sum/3,
4 matrix_diff/3,
5 matrix_mult_scal/3,
6 matrix_div_scal/3,
7 dot_product/3,
8 cholesky_decomposition/2,
9 matrix_inversion/2,
10 matrix_inv_triang/2,
11 matrix_identity/2,
12 matrix_diagonal/2,
13 determinant/2,
14 list0/2
15 ]).
39:- use_module(library(clpfd), [transpose/2]).
43matrix_div_scal(A,V,B):- 44 maplist(maplist(div(V)),A,B). 45 46div(A,B,C):- 47 C is B/A.
51matrix_mult_scal(A,V,B):- 52 maplist(maplist(mult(V)),A,B). 53 54mult(A,B,C):- 55 C is A*B.
?- determinant([[2,-1,0],[-1,2,-1],[0,-1,2]],D). D = 3.999999999999999.
63determinant(A,Det):- 64 cholesky_decomposition(A,L), 65 get_diagonal(L,D), 66 foldl(prod,D,1,DetL), 67 Det is DetL*DetL. 68 69prod(A,P0,P):- 70 P is P0*A.
76matrix_identity(N,M):- 77 numlist(1,N,L), 78 maplist(row(N,L),L,M). 79 80 81row(N,L,I,Row):- 82 length(Row,N), 83 maplist(elem(I),L,Row). 84 85elem(I,IP,El):- 86 (I==IP-> 87 El=1 88 ; 89 El=0 90 ).
diag(Vect)
97matrix_diagonal(Vect,M):- 98 length(Vect,N), 99 numlist(1,N,L), 100 maplist(row_diag(N,L),L,Vect,M). 101 102 103row_diag(N,L,I,El,Row):- 104 length(Row,N), 105 maplist(elem_diag(I,El),L,Row). 106 107elem_diag(I,Comp,IP,El):- 108 (I==IP-> 109 El=Comp 110 ; 111 El=0 112 ).
?- matrix_inversion([[2,-1,0],[-1,2,-1],[0,-1,2]],L). L = [[0.7499999999999999, 0.5000000000000001, 0.2500000000000001], [0.5000000000000001, 1.0000000000000004, 0.5000000000000002], [0.2500000000000001, 0.5000000000000002, 0.7500000000000001]].
123matrix_inversion(A,B):-
124 cholesky_decomposition(A,L),
125 matrix_inv_triang(L,LI),
126 transpose(LI,LIT),
127 matrix_multiply(LIT,LI,B).
?- matrix_inv_triang([[2,0,0],[-1,2,0],[0,-1,2]],L). L = [[0.5, 0.0, 0.0], [0.25, 0.5, 0.0], [0.125, 0.25, 0.5]].
139matrix_inv_triang(L1,L2):- 140 get_diagonal(L1,D), 141 maplist(inv,D,ID), 142 diag(ID,IDM), 143 matrix_multiply(IDM,L1,LL1), 144 list_to_term(LL1,LT), 145 length(LL1,N), 146 matrix_inv_i(1,N,LT), 147 term_to_list(LT,N,LL2), 148 matrix_multiply(LL2,IDM,L2). 149 150 151matrix_inv_i(N,N,_LT):-!. 152 153matrix_inv_i(I,N,LT):- 154 matrix_inv_j(0,I,N,LT), 155 I1 is I+1, 156 matrix_inv_i(I1,N,LT). 157 158matrix_inv_j(I,I,_N,_LT):-!. 159 160matrix_inv_j(J,I,N,LT):- 161 get_v(I,J,N,LT,Vij), 162 V_ij is -Vij, 163 set_v(I,J,N,LT,V_ij), 164 J1 is J+1, 165 matrix_inv_k(J1,J,I,N,LT), 166 matrix_inv_j(J1,I,N,LT). 167 168matrix_inv_k(I,_J,I,_N,_LT):-!. 169 170matrix_inv_k(K,J,I,N,LT):- 171 get_v(I,K,N,LT,Vik), 172 get_v(K,J,N,LT,Vkj), 173 get_v(I,J,N,LT,Vij), 174 NVij is Vij-Vik*Vkj, 175 set_v(I,J,N,LT,NVij), 176 K1 is K+1, 177 matrix_inv_k(K1,J,I,N,LT). 178 179 180diag(L,D):- 181 length(L,N), 182 NN is N*N, 183 list0(NN,L0), 184 DT =..[a|L0], 185 N1 is N-1, 186 numlist(0,N1,In), 187 maplist(set_diag(DT,N),In,L), 188 term_to_list(DT,N,D). 189 190set_diag(D,N,I,V):- 191 set_v(I,I,N,D,V). 192 193inv(A,B):- 194 B is 1.0/A. 195 196get_diagonal(L,D):- 197 length(L,N), 198 list_to_term(L,LT), 199 get_diag(0,N,LT,D). 200 201get_diag(N,N,_L,[]):-!. 202 203get_diag(N0,N,L,[H|R]):- 204 get_v(N0,N0,N,L,H), 205 N1 is N0+1, 206 get_diag(N1,N,L,R). 207 208list_to_term(L,LT):- 209 append(L,LL), 210 LT=..[a|LL].
?- matrix_multiply([[1,2],[3,4],[5,6]], [[1,1,1],[1,1,1]],R). R = [[3, 3, 3], [7, 7, 7], [11, 11, 11]].
code from http://stackoverflow.com/questions/34206275/matrix-multiplication-with-prolog
219matrix_multiply(X,Y,M) :- 220 matrix_mul(X,Y,M0), 221 maplist(maplist(is),M,M0). 222 223matrix_mul(X,Y,M) :- 224 transpose(Y,T), 225 maplist(row_multiply(T),X,M). 226 227row_multiply(T,X,M) :- 228 maplist(dot_product(X),T,M).
233dot_product([X|Xs],[T|Ts],M) :- 234 foldl(mul,Xs,Ts,X*T,M). 235mul(X,T,M,M+X*T).
238matrix_diff(X,Y,S):- 239 maplist(maplist(diff),X,Y,S). 240 241diff(A,B,C):- 242 C is A-B.
matrix_sum([[1,2],[3,4],[5,6]],[[1,2],[3,4],[5,6]],M). M = [[2, 4], [6, 8], [10, 12]].
248matrix_sum(X,Y,S):- 249 maplist(maplist(sum),X,Y,S). 250 251sum(A,B,C):- 252 C is A+B.
cholesky_decomposition([[25, 15, -5], [15, 18, 0], [-5, 0, 11]],L). L = [[5.0, 0, 0], [3.0, 3.0, 0], [-1.0, 1.0, 3.0]]. cholesky_decomposition([[18, 22, 54, 42],[22, 70, 86, 62],[ 54, 86, 174, 134],[ 42, 62, 134, 106]],L). L = [[4.242640687119285, 0, 0, 0], [5.185449728701349, 6.565905201197403, 0, 0], [12.727922061357857, 3.0460384954008553, 1.6497422479090704, 0], [9.899494936611667, 1.624553864213788, 1.8497110052313648, 1.3926212476456026]].
263cholesky_decomposition(A,L):- 264 append(A,AL), 265 length(AL,NL), 266 AM =..[a|AL], 267 list0(NL,LL), 268 LM=..[l|LL], 269 length(A,N), 270 cholesky_i(0,N,AM,LM), 271 term_to_list(LM,N,L). 272 273cholesky_i(N,N,_A,_L):-!. 274 275cholesky_i(I,N,A,L):- 276 cholesky_j(0,I,N,A,L), 277 I1 is I+1, 278 cholesky_i(I1,N,A,L). 279 280cholesky_j(I,I,N,A,L):-!, 281 cholesky_k(0,I,I,N,0,S,L), 282 get_v(I,I,N,A,Aii), 283 V is sqrt(Aii-S), 284 set_v(I,I,N,L,V). 285 286 287cholesky_j(J,I,N,A,L):- 288 cholesky_k(0,J,I,N,0,S,L), 289 get_v(I,J,N,A,Aij), 290 get_v(J,J,N,L,Ljj), 291 V is 1.0/Ljj*(Aij-S), 292 set_v(I,J,N,L,V), 293 J1 is J+1, 294 cholesky_j(J1,I,N,A,L). 295 296cholesky_k(J,J,_I,_N,S,S,_L):-!. 297 298cholesky_k(K,J,I,N,S0,S,L):- 299 get_v(I,K,N,L,Lik), 300 get_v(J,K,N,L,Ljk), 301 S1 is S0+Lik*Ljk, 302 K1 is K+1, 303 cholesky_k(K1,J,I,N,S1,S,L). 304 305get_v(I,J,N,M,V):- 306 Argij is I*N+J+1, 307 arg(Argij,M,V). 308 309set_v(I,J,N,M,V):- 310 Argij is I*N+J+1, 311 setarg(Argij,M,V). 312 313 314term_to_list(T,N,L):- 315 T=..[_|E], 316 identify_rows(E,N,L). 317 318identify_rows([],_N,[]):-!. 319 320identify_rows(E,N,[R|L]):- 321 length(R,N), 322 append(R,Rest,E), 323 identify_rows(Rest,N,L).
327list0(0,[]):-!. 328 329list0(N,[0|T]):- 330 N1 is N-1, 331 list0(N1,T)
matrix
This module performs matrix operations. Impemented operations:
The library was developed for dealing with multivariate Gaussian distributions, that's the reson for the focus on positive semi-definite matrices