Matrix Operations

Matrices are used for solving linear systems of equations and are used extensively in polynomial fitting, adaptive equalization, and filter design. In liquid , matrices are represented as just arrays of a single dimension, and do not rely on special objects for their manipulation. This is to help portability of the code and ease of integration into other libraries. Here is a simple example of the matrix interface:


#include <liquid/liquid.h>

int main() {
    // designate X as a 4 x 4 matrix
    float X[16] = {
       0.84382,  -2.38304,   1.43061,  -1.66604,
       3.99475,   0.88066,   4.69373,   0.44563,
       7.28072,  -2.06608,   0.67074,   9.80657,
       6.07741,  -3.93099,   1.22826,  -0.42142};
    matrixf_print(X,4,4);

    // L/U decomp (Doolittle's method)
    float L[16], U[16], P[16];
    matrixf_ludecomp_doolittle(X,4,4,L,U,P);
}

Notice that all routines for the type float are prefaced with matrixf . This follows the naming convention of the standard C library routines which append an f to the end of methods operating on floating-point precision types. Similar matrix interfaces exist in liquid for double ( matrix ), double complex ( matrixc ), and float complex ( matrixcf ).

Basic math operations

This section describes the basic matrix math operations, including addition, subtraction, point-wise multiplication and division, transposition, and initializing the identity matrix.

matrix_access (access element)

Because matrices in liquid are really just one-dimensional arrays, indexing matrix values for storage or retrieval is as straightforward as indexing the array itself. liquid also provides a simple macro for ensuring the proper value is returned. matrix_access(X,R,C,r,c) will access the element of a \(R \times C\) matrix \(\vec{X}\) at row \(r\) and column \(c\) . This method is really just a pre-processor macro which performs a literal string replacement


#define matrix_access(X,R,C,r,c) ((X)[(r)*(C)+(c)])

and can be used for both setting and retrieving values of a matrix. For example,


X =
  0.406911   0.118444   0.923281   0.827254   0.463265
  0.038897   0.132381   0.061137   0.880045   0.570341
  0.151206   0.439508   0.695207   0.215935   0.999683
  0.808384   0.601597   0.149171   0.975722   0.205819

float v = matrix_access(X,4,5,0,1);
v =
  0.118444

matrix_access(X,4,5,2,3) = 0;
X =
  0.406911   0.118444   0.923281   0.827254   0.463265
  0.038897   0.132381   0.061137   0.880045   0.570341
  0.151206   0.439508   0.695207   0.0        0.999683
  0.808384   0.601597   0.149171   0.975722   0.205819

Because this method is really just a macro, there is no error-checking to ensure that one is accessing the matrix within its memory bounds. Therefore, special care must be taken when programming. Furthermore, matrix_access() can be used for all matrix types ( matrixf , matrixcf , etc.).

matrixf_add, matrixf_sub, matrixf_pmul, and matrixf_pdiv (scalar math operations)

The matrixf_add(*x,*y,*z,m,n) , matrixf_sub(*x,*y,*z,m,n) , matrixf_pmul(*x,*y,*z,m,n) , and matrixf_pdiv(*x,*y,*z,m,n) methods perform point-wise (scalar) addition, subtraction, multiplication, and division of the elements of two \(n \times m\) matrices, \(\vec{X}\) and \(\vec{Y}\) . That is, \(\vec{Z}_{i,k} = \vec{X}_{i,k} + \vec{Y}_{i,k}\) for all \(i\) , \(k\) . The same holds true for subtraction, multiplication, and division. It is very important to understand the difference between the methods matrixf_pmul() and matrixf_mul() , as well as matrixf_pdiv() and matrixf_div() . In each case the latter performs a vastly different operation from matrixf_mul() and matrixf_div() (see Sections [section-matrix-mul] and [section-matrix-div] , respectively).


X =                         Y =
  0.59027   0.83429           0.764108   0.741641
  0.67779   0.19793           0.660932   0.041723
  0.95075   0.33980           0.972282   0.347090

matrixf_pmul(X,Y,Z,2,3);
Z =
  0.4510300   0.6187437
  0.4479731   0.0082582
  0.9243971   0.1179412

matrixf_trans(), matrixf_hermitian() (transpose matrix)

The matrixf_trans(X,m,n,XT) method performs the conjugate matrix transpose operation on an \(m \times n\) matrix \(\vec{X}\) . That is, the matrix is flipped on its main diagonal and the conjugate of each element is taken. Formally, \(\vec{A}^T_{i,j} = \vec{A}_{j,i}^*\) . Here's a simple example:

$$ \left[ \begin{array}{ccc} 0 & 1 & 2 \\ 3 & 4 & 5 \\ \end{array} \right]^T = \left[ \begin{array}{cc} 0 & 3 \\ 1 & 4 \\ 2 & 5 \\ \end{array} \right] $$

Similarly, the matrixf_hermitian(X,m,n,XH) computes the Hermitian transpose which is identical to the regular transpose but without the conjugation operation, viz \(\vec{A}^H_{i,j} = \vec{A}_{j,i}\) .

matrixf_eye() (identity matrix)

The matrixf_eye(*x,n) method generates the \(n \times n\) identity matrix \(\vec{I}_n\) :

$$ \vec{I}_n = \begin{bmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 \\ \\ 0 & 0 & \cdots & 1 \\ \end{bmatrix} $$

Elementary math operations

This section describes elementary math operations for linear systems of equations.

matrixf_swaprows() (swap rows)

Matrix row-swapping is often necessary to express a matrix in its row-reduced echelon form. The matrixf_swaprows(*X,m,n,i,j) method simply swaps rows \(i\) and \(j\) of an \(m \times n\) matrix \(\vec{X}\) , viz


x =
  0.84381998 -2.38303995  1.43060994 -1.66603994
  3.99475002  0.88066000  4.69372988  0.44563001
  7.28072023 -2.06608009  0.67074001  9.80657005
  6.07741022 -3.93098998  1.22826004 -0.42142001

matrixf_swaprows(x,4,4,0,2);
  7.28072023 -2.06608009  0.67074001  9.80657005
  3.99475002  0.88066000  4.69372988  0.44563001
  0.84381998 -2.38303995  1.43060994 -1.66603994
  6.07741022 -3.93098998  1.22826004 -0.42142001

matrixf_pivot() (pivoting)

NOTE--terminology for "pivot" is different from literature.

Given an \(n \times m\) matrix \(\vec{A}\) ,

$$ \vec{A} = \begin{bmatrix} A_{0,0} & A_{0,1} & \cdots & A_{0,m-1} \\ A_{1,0} & A_{1,1} & \cdots & A_{1,m-1} \\ \\ A_{n-1,0} & A_{n-1,1} & \cdots & A_{n-1,m-1} \end{bmatrix} $$

pivoting \(\vec{A}\) around \(\vec{A}_{a,b}\) gives

$$ \vec{B}_{i,j} = \left( \frac{\vec{A}_{i,b}}{\vec{A}_{a,b}} \right) \vec{A}_{a,j} - \vec{A}_{i,j} \forall i \ne a $$

The pivot element must not be zero. Row \(a\) is left unchanged in \(\vec{B}\) . All elements of \(\vec{B}\) in column \(b\) are zero except for row \(a\) . This is accomplished in liquid with the matrixf_pivot(*A,m,n,i,j) method. For our example \(4 \times 4\) matrix \(\vec{x}\) , pivoting around \(\vec{x}_{1,2}\) gives:


matrixf_pivot(x,4,4,1,2);
  0.37374675  2.65145779  0.00000000  1.80186427
  3.99475002  0.88066000  4.69372988  0.44563001
 -6.70986557  2.19192743  0.00000000 -9.74288940
 -5.03205967  4.16144180  0.00000000  0.53803295

matrixf_mul() (multiplication)

Multiplication of two input matrices \(\vec{A}\) and \(\vec{B}\) is accomplished with the matrixf_mul(*A,ma,na,*B,mb,nb,*C,mc,nc) method, and is not to be confused with matrixf_pmul() in [section-matrix-mathop] . If \(\vec{A}\) is \(m \times n\) and \(\vec{B}\) is \(n \times p\) , then their product is computed as

$$ \bigl( \vec{A} \vec{B} \bigr)_{i,j} = \sum_{r=0}^{n-1} { \vec{A}_{i,r} \vec{B}_{r,j} } $$

Note that the number of columns of \(\vec{A}\) must be equal to the number of rows of \(\vec{B}\) , and that the resulting matrix is of size \(m \times p\) (the number of rows in \(\vec{A}\) and columns in \(\vec{B}\) ).


A =                 B =
    1   2   3           1   2   3
    4   5   6           4   5   6
                        7   8   9
matrixf_mul(A,2,3, B,3,3, C,2,3);

C =
    30  36  42
    66  81  96

Transpose multiplication

liquid also implements transpose-multiplication operations on an \(m \times n\) matrix \(\vec{X}\) , commonly used in signal processing. [section-matrix-trans] describes the difference between the \(\left(\cdot\right)^T\) and \(\left(\cdot\right)^H\) operations. The interface for transpose-multiplications in liquid is tabulated below for an input \(m \times n\) matrix \(\vec{X}\) .


.. table
_operation_             & _output dimensions_   & _interface_\\\otoprule
$\vec{X}  \vec{X}^T$    & $m \times m$          & `matrixcf_mul_transpose(x,m,n,xxT)`
$\vec{X}  \vec{X}^H$    & $m \times m$          & `matrixcf_mul_hermitian(x,m,n,xxH)`
$\vec{X}^T\vec{X}  $    & $n \times n$          & `matrixcf_transpose_mul(x,m,n,xTx)`
$\vec{X}^H\vec{X}  $    & $n \times n$          & `matrixcf_transpose_mul(x,m,n,xHx)`

Complex math operations

More complex math operations are described here, including matrix inversion, square matrix determinant, Gauss-Jordan elimination, and lower/upper decomposition routines using both Crout's and Doolittle's methods.

matrixf_inv() (inverse)

Matrix inversion is accomplished with the matrixf_inv(*X,m,n) method.


.. footnote
While matrix inversion requires a square matrix, `liquid`
internally checks to ensure $m=n$ on the input size for
$\vec{X}$.

Given an \(n \times n\) matrix \(\vec{A}\) , liquid augments with \(\vec{I}_n\) :

$$ \left[\vec{A}|\vec{I}_n\right] = \left[ \begin{array}{cccc|cccc} A_{0,0} & A_{0,1} & \cdots & A_{0,m-1} & 1 & 0 & \cdots & 0 \\ A_{1,0} & A_{1,1} & \cdots & A_{1,m-1} & 0 & 1 & \cdots & 0 \\ & & & & & & & \\ A_{n-1,0} & A_{n-1,1} & \cdots & A_{n-1,m-1} & 0 & 0 & \cdots & 1 \\ \end{array} \right] $$

Next liquid performs elementary operations to convert to its row-reduced echelon form. The resulting matrix has the identity matrix on the left and \(\vec{A}^{-1}\) on its right, viz

$$ \left[\vec{I}_n|\vec{A}^{-1}\right] = \left[ \begin{array}{cccc|cccc} 1 & 0 & \cdots & 0 & A^{-1}_{0,0} & A^{-1}_{0,1} & \cdots & A^{-1}_{0,m-1} \\ 0 & 1 & \cdots & 0 & A^{-1}_{1,0} & A^{-1}_{1,1} & \cdots & A^{-1}_{1,m-1} \\ & & & & & & & \\ 0 & 0 & \cdots & 1 & A^{-1}_{n-1,0} & A^{-1}_{n-1,1} & \cdots & A^{-1}_{n-1,m-1} \\ \end{array} \right] $$

The matrixf_inv() method uses Gauss-Jordan elimination (see matrixf_gjelim() ) for row reduction and back-substitution. Pivot elements in \(\vec{A}\) with the largest magnitude are chosen to help stability in floating-point arithmetic.


matrixf_inv(x,4,4);
 -0.33453920  0.04643385 -0.04868321  0.23879384
 -0.42204019  0.12152659 -0.07431178  0.06774280
  0.35104612  0.15256262  0.04403552 -0.20177667
  0.13544561 -0.01930523  0.11944833 -0.14921521

matrixf_div()

The matrixf_div(*X,*Y,*Z,*n) method simply computes \(\vec{Z} = \vec{Y}^{-1}\vec{X}\) where \(\vec{X}\) , \(\vec{Y}\) , and \(\vec{Z}\) are all \(n \times n\) matrices.

matrixf_linsolve() (solve linear system of equations)

The matrixf_linsolve(*A,n,*b,*x,opts) method solves a set of \(n\) linear equations \(A\vec{x} = \vec{b}\) where \(A\) is an \(n \times n\) matrix, and \(\vec{x}\) and \(\vec{b}\) are \(n \times 1\) vectors. The opts argument is reserved for future development and can be ignored by setting to NULL .

matrixf_cgsolve() (solve linear system of equations)

The matrixf_cgsolve(*A,n,*b,*x,opts) method solves \(Ax = b\) using the conjugate gradient method where \(A\) is an \(n \times n\) symmetric positive-definite matrix. The opts argument is reserved for future development and can be ignored by setting to NULL . Listed below is a basic example:


A =
  2.9002075   0.1722705   1.3046706   1.8082311
  0.1722705   1.0730995   0.2497573   0.1470398
  1.3046706   0.2497573   0.8930279   1.1471686
  1.8082311   0.1470398   1.1471686   1.5155975
b =
 11.7622252
 -1.0541668
  5.7372437
  8.1291904
matrixf_cgsolve(A,4,4, x_hat, NULL)
x_hat =
  2.8664699
 -1.8786657
  1.1224079
  1.2764599

For a more complete example, see examples/cgsolve_example.c located under the main project directory.

matrixf_det() (determinant)

The matrixf_det(*X,m,n) method computes the determinant of an \(n \times n\) matrix \(\vec{X}\) . In liquid , the determinant is computed by L/U decomposition of \(\vec{A}\) using Doolittle's method (see matrixf_ludecomp_doolittle ) and then computing the product of the diagonal elements of \(\vec{U}\) , viz

$$ \det\left(\vec{A}\right) = \left|\vec{A}\right| = \prod_{k=0}^{n-1}{\vec{U}_{k,k}} $$

This is equivalent to performing L/U decomposition using Crout's method and then computing the product of the diagonal elements of \(\vec{L}\) .


matrixf_det(X,4,4) = 585.40289307

matrixf_ludecomp_crout() (LU Decomposition, Crout's Method)

Crout's method decomposes a non-singular \(n\times n\) matrix \(\vec{A}\) into a product of a lower triangular \(n \times n\) matrix \(\vec{L}\) and an upper triangular \(n \times n\) matrix \(\vec{U}\) . In fact, \(\vec{U}\) is a unit upper triangular matrix (its values along the diagonal are 1). The matrixf_ludecomp_crout(*A,m,n,*L,*U,*P) implements Crout's method.

$$ \vec{L}_{i,k} = \vec{A}_{i,k} - \sum_{t=0}^{k-1}{ \vec{L}_{i,t} \vec{U}_{t,k} } \ \forall k \in \{0,n-1\}, i \in \{k,n-1\} $$ $$ \vec{U}_{k,j} = \left[ \vec{A}_{k,j} - \sum_{t=0}^{k-1}{ \vec{L}_{k,t} \vec{U}_{t,j} } \right] / \vec{L}_{k,k} \ \forall k \in \{0,n-1\}, j \in \{k+1,n-1\} $$

matrixf_ludecomp_crout(X,4,4,L,U,P)
L =
  0.84381998  0.00000000  0.00000000  0.00000000
  3.99475002 12.16227055  0.00000000  0.00000000
  7.28072023 18.49547005 -8.51144791  0.00000000
  6.07741022 13.23228073 -6.81350422 -6.70173073
U =
  1.00000000 -2.82410932  1.69539714 -1.97440207
  0.00000000  1.00000000 -0.17093502  0.68514121
  0.00000000  0.00000000  1.00000000 -1.35225296
  0.00000000  0.00000000  0.00000000  1.00000000

matrixf_ludecomp_doolittle() (LU Decomposition, Doolittle's Method)

Doolittle's method is similar to Crout's except it is the lower triangular matrix that is left with ones on the diagonal. The update algorithm is similar to Crout's but with a slight variation: the upper triangular matrix is computed first. The matrixf_ludecomp_doolittle(*A,m,n,*L,*U,*P) implements Doolittle's method.

$$ \vec{U}_{k,j} = \vec{A}_{k,j} - \sum_{t=0}^{k-1}{ \vec{L}_{k,t} \vec{U}_{t,j} } \ \forall k \in \{0,n-1\}, j \in \{k,n-1\} $$ $$ \vec{L}_{i,k} = \left[ \vec{A}_{i,k} - \sum_{t=0}^{k-1}{ \vec{L}_{i,t} \vec{U}_{t,k} } \right] / \vec{U}_{k,k} \ \forall k \in \{0,n-1\}, i \in \{k+1,n-1\} $$

Here is a simple example:


matrixf_ludecomp_doolittle(X,4,4,L,U,P)
L =
  1.00000000  0.00000000  0.00000000  0.00000000
  4.73412609  1.00000000  0.00000000  0.00000000
  8.62828636  1.52072513  1.00000000  0.00000000
  7.20225906  1.08797777  0.80051047  1.00000000
U =
  0.84381998 -2.38303995  1.43060994 -1.66603994
  0.00000000 12.16227150 -2.07895803  8.33287334
  0.00000000  0.00000000 -8.51144791 11.50963116
  0.00000000  0.00000000  0.00000000 -6.70172977

matrixf_qrdecomp_gramschmidt() (QR Decomposition, Gram-Schmidt algorithm)

liquid implements Q/R decomposition with the matrixf_qrdecomp_gramschmidt(*A,m,n,*Q,*R) method which factors a non-singular \(n \times n\) matrix \(\vec{A}\) into product of an orthogonal matrix \(\vec{Q}\) and an upper triangular matrix \(\vec{R}\) , each \(n \times n\) . That is, \(\vec{A} = \vec{Q}\vec{R}\) where \(\vec{Q}^T\vec{Q} = \vec{I}_n\) and \(\vec{R}_{i,j} = 0\,\, \forall_{i > j}\) . Building on the previous example for our test \(4 \times 4\) matrix \(\vec{X}\) , the Q/R factorization is


matrixf_qrdecomp_gramschmidt(X,4,4,Q,R)
Q =
  0.08172275 -0.57793844  0.57207584  0.57622749
  0.38688579  0.63226062  0.66619849 -0.08213031
  0.70512730  0.13563085 -0.47556636  0.50816941
  0.58858842 -0.49783322  0.05239720 -0.63480729
R =
 10.32539940 -3.62461853  3.12874746  6.70309162
  0.00000000  3.61081028  1.62036073  2.78449297
  0.00000000  0.00000000  3.69074893 -5.34197950
  0.00000000  0.00000000  0.00000000  4.25430155

matrixf_chol() (Cholesky Decomposition)

Compute Cholesky decomposition of an \(n \times n\) symmetric/Hermitian positive-definite matrix as \(\vec{A} = \vec{L}\vec{L}^T\) where \(\vec{L}\) is \(n \times n\) and lower triangular. An \(n \times n\) matrix is positive definite if \(\Re\{v^T \vec{A} v\} > 0\) for all non-zero vectors \(v\) . Note that \(\vec{A}\) can be either complex or real. Shown below is an example of the Cholesky decomposition of a \(4 \times 4\) positive definite real matrix.


A =
  1.0201000  -1.4341999   0.3232000  -1.0302000
 -1.4341999   2.2663999   0.5506001   1.2883999
  0.3232000   0.5506001   4.2325001  -1.4646000
 -1.0302000   1.2883999  -1.4646000   5.0101995
matrixf_chol(A,4,Lp)
  1.0100000   0.0000000   0.0000000   0.0000000
 -1.4200000   0.5000000   0.0000000   0.0000000
  0.3200000   2.0100000   0.3000003   0.0000000
 -1.0200000  -0.3199999  -1.6499993   1.0700010

matrixf_gjelim() (Gauss-Jordan Elimination)

The matrixf_gjelim(*X,m,n) method in liquid performs the Gauss-Jordan elimination on a matrix \(\vec{X}\) . Gauss-Jordan elimination converts a \(m \times n\) matrix into its row-reduced echelon form using elementary matrix operations (e.g. pivoting). This can be used to solve a linear system of \(n\) equations \(\vec{A}\vec{x} = \vec{b}\) for the unknown vector \(\vec{x}\)

$$ \begin{bmatrix} A_{0,0} & A_{0,1} & \cdots & A_{0,n-1} \\ A_{1,0} & A_{1,1} & \cdots & A_{1,n-1} \\ \\ A_{n-1,0} & A_{n-1,1} & \cdots & A_{n-1,n-1} \end{bmatrix} \begin{bmatrix} x_{0} \\ x_{1} \\ \\ x_{n-1} \end{bmatrix} = \begin{bmatrix} b_{0} \\ b_{1} \\ \\ b_{n-1} \end{bmatrix} $$

The solution for \(\vec{x}\) is given by inverting \(\vec{A}\) and multiplying by \(\vec{b}\) , viz

$$ \vec{x} = \vec{A}^{-1}\vec{b} $$

This is also equivalent to augmenting \(\vec{A}\) with \(\vec{b}\) and converting it to its row-reduced echelon form. If \(\vec{A}\) is non-singular the resulting \(n \times n+1\) matrix will hold \(\vec{x}\) in its last column. The row-reduced echelon form of a matrix is computed in liquid using the Gauss-Jordan elimination algorithm, and can be invoked as such:


Ab =
  0.84381998 -2.38303995  1.43060994 -1.66603994  0.91488999
  3.99475002  0.88066000  4.69372988  0.44563001  0.71789002
  7.28072023 -2.06608009  0.67074001  9.80657005  1.06552994
  6.07741022 -3.93098998  1.22826004 -0.42142001 -0.81707001
matrixf_gjelim(Ab,4,5)
  1.00000000 -0.00000000  0.00000000 -0.00000000 -0.51971692
 -0.00000000  1.00000000  0.00000000  0.00000000 -0.43340963
 -0.00000000 -0.00000000  1.00000000 -0.00000000  0.64247853
  0.00000000 -0.00000000 -0.00000000  0.99999994  0.35925382

Notice that the result contains \(\vec{I}_n\) in its first \(n\) rows and \(n\) columns (to within machine precision).


.. footnote
row permutations (swapping) might have occurred.