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    ]).

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

author
- Fabrizio Riguzzi
license
- Artistic License 2.0
   39:- use_module(library(clpfd), [transpose/2]).
 matrix_div_scal(+A, +V, -B) is det
divide matrix A by scalar V
   43matrix_div_scal(A,V,B):-
   44  maplist(maplist(div(V)),A,B).
   45
   46div(A,B,C):-
   47  C is B/A.
 matrix_mult_scal(+A, +V, -B) is det
multiply matrix A by scalar V
   51matrix_mult_scal(A,V,B):-
   52  maplist(maplist(mult(V)),A,B).
   53
   54mult(A,B,C):-
   55  C is A*B.
 determinant(+A, -D) is det
computes the determinant for a positive semi-definite matrix. Uses the Cholenski decomposition
?- 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.
 matrix_identity(+N, -M) is det
generates identity matrix of size NxN
   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  ).
 matrix_diagonal(+Vect, -M) is det
generates diagonal matrix from vector Vect M=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(+M, -IM) is det
inversion of a positive semi-definite matrix. Uses the Cholenski decomposition
?- 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(+M, -IM) is det
inversion of a lower triangular matrix code from http://www.mymathlib.com/c_source/matrices/linearsystems/unit_lower_triangular.c http://www.mcs.csueastbay.edu/~malek/TeX/Triangle.pdf code from
?- 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(+X, +Y, -M) is det
X(N*P),Y(P*M),M(N*M)
?- 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).
 dot_product(+X, +Y, -D) is det
computes the dot produce of two vectors
  233dot_product([X|Xs],[T|Ts],M) :-
  234    foldl(mul,Xs,Ts,X*T,M).
  235mul(X,T,M,M+X*T).
 matrix_diff(+A, +B, -C) is det
  238matrix_diff(X,Y,S):-
  239  maplist(maplist(diff),X,Y,S).
  240
  241diff(A,B,C):-
  242  C is A-B.
 matrix_sum(+A, +B, -C) is det
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(+A, -L) is det
computes the Cholesky decomposition of a positive semi-definite matrix code from https://rosettacode.org/wiki/Cholesky_decomposition#C
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).
 list0(+N, -L) is det
returns a list of N zeros
  327list0(0,[]):-!.
  328
  329list0(N,[0|T]):-
  330  N1 is N-1,
  331  list0(N1,T)