# cl-simple-matrix

**Repository Path**: repl-top/cl-simple-matrix

## Basic Information

- **Project Name**: cl-simple-matrix
- **Description**: No description available
- **Primary Language**: Unknown
- **License**: GPL-3.0
- **Default Branch**: main
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 0
- **Forks**: 0
- **Created**: 2025-07-30
- **Last Updated**: 2025-07-30

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

#+TITLE: simple-matrix
#+AUTHOR: Guillaume LE VAILLANT
#+DATE: 2025-06-15
#+EMAIL: glv@posteo.net
#+LANGUAGE: en
#+OPTIONS: num:nil toc:nil html-postamble:nil html-scripts:nil
#+HTML_DOCTYPE: html5

* Description
*simple-matrix* is a Common Lisp library implementing some functions to
work with matrices.

This library is designed to be generic, to be used for a few simple
computations. If you need something to work very fast on a lot of big matrices,
this is not the library you are looking for.

Matrices are represented as arrays with two dimensions containing elements of
type ~number~. So you can have a matrix containing integers, rationals, floats
and complexes at the same time if you want.

* License
*simple-matrix* is released under the GPL-3 license or later. See the
[[file:LICENSE][LICENSE]] file for details.

* API
After putting the library in a place where your Common Lisp system can see it,
it can be loaded with the usual:
#+BEGIN_SRC lisp
  (asdf:load-system "simple-matrix")
#+END_SRC

As the code is in a single file, you could also just load it directly:
#+BEGIN_SRC lisp
  (load "/path/to/simple-matrix.lisp")
#+END_SRC

The functions will then be available in the ~simple-matrix~ package.

** Parameters
Numbers with an amplitude lesser than the ~*epsilon*~ parameter will be
considered equal to 0 by some functions (1.0d-6 by default).

The ~*qr-variant*~ parameter indicates which algorithm to use for the QR
decomposition:
 - ~:householder~ to use Householder transformations and return a full QR
   decomposition. This is the variant used by default.
 - ~:householder-reduced~ to use Householder transformations and return
   a reduced QR decomposition.
 - ~:givens~ to use Givens rotations and return a full QR decomposition.
   This can make computations faster for matrices containing a lot of zeros.
 - ~:givens-reduced~ to use Givens rotations and return a reduced QR
   decomposition.
 - ~:gram-schmidt~ to use the Gram-Schmidt process and return a reduced QR
   decomposition.
The reduced QR variants can be useful, for example, to make computations faster
for the least squares method on a matrix with a lot more rows than columns, or
to compute a reduced singular value decomposition.

The ~*schur-variant*~ parameter indicates which variant of the Schur
decomposition to use. It can be ~:complex~ (the default), or ~:real~.

The ~*max-iterations*~ parameter indicates the maximum number of iterations
used to compute some functions (like ~mexp~, ~msqrt~, ~mschur~, ~meigen~, ~mqz~
or ~msvd~). It is 100 by default. It can be ~nil~ to make the computation stop
only when the function considers it has converged (which could be never).

The ~*diagonalize-p*~ parameter indicates if some functions (like ~mexp~ or
~msqrt~) will try to diagonalize the matrix to compute their result. It is ~t~
by default.

** Functions
*** Matrix creation and element access
As matrices are represented as arrays with two dimensions, their elements can
be accessed using the regular ~aref~ function.

#+BEGIN_SRC lisp
  (matrix r c &optional elements by-column) => matrix
#+END_SRC
Create a matrix with /r/ rows and /c/ columns. If the /elements/ are provided,
use them to fill the matrix either row by row if /by-column/ is ~nil~ (the
default), or column by column if /by-column/ is ~t~. The /elements/ can be
provided as a list, a vector, a row matrix or a column matrix. It the
/elements/ are not provided, the matrix will be full of zeros.

#+BEGIN_SRC lisp
  (vec->mat v) => matrix
#+END_SRC
Convert the vector /v/ to a column matrix.

#+BEGIN_SRC lisp
  (mat->vec m) => vector
#+END_SRC
Convert the column matrix /M/ to a vector.

#+BEGIN_SRC lisp
  (mvec m) => matrix
#+END_SRC
Return a column matrix containing the concatenated columns of the matrix /M/.

#+BEGIN_SRC lisp
  (mi n) => matrix
#+END_SRC
Create an identity matrix with /n/ rows and /n/ columns.

#+BEGIN_SRC lisp
  (mcopy m) => matrix
#+END_SRC
Return a copy of the matrix /M/.

#+BEGIN_SRC lisp
  (mround m) => matrix
#+END_SRC
Return a copy of the matrix /M/ with its elements rounded to the precision
defined by ~*epsilon*~.

#+BEGIN_SRC lisp
  (msub m &key row-start row-end column-start column-end) => matrix
#+END_SRC
Return a submatrix of /M/, starting at row /row-start/ and column
/column-start/, and ending before row /row-end/ and column /column-end/.
This function is setf-able.

#+BEGIN_SRC lisp
  (mrow m i) => vector
#+END_SRC
Return row /i/ of matrix /M/ as a vector.
This function is setf-able.

#+BEGIN_SRC lisp
  (mcol m j) => vector
#+END_SRC
Return column /j/ of matrix /M/ as a vector.
This function is setf-able.

#+BEGIN_SRC lisp
  (mdiag m) => vector
#+END_SRC
Return diagonal of matrix /M/ as a vector.
This function is setf-able.

#+BEGIN_SRC lisp
  (mminor m i j) => matrix
#+END_SRC
Return the matrix /M/ with row /i/ and column /j/ removed.

*** Predicates
These functions take ~*epsilon*~ into consideration.

#+BEGIN_SRC lisp
  (m= m &rest matrices) => boolean
#+END_SRC
Return ~t~ if the matrices are equal, and ~nil~ otherwise.

#+BEGIN_SRC lisp
  (mdiagonal-p m) => boolean
#+END_SRC
Return ~t~ if /M/ is a diagonal matrix, and ~nil~ otherwise.

#+BEGIN_SRC lisp
  (mupper-triangular-p m) => boolean
#+END_SRC
Return ~t~ if the matrix /M/ is in upper triangular form, and ~nil~ otherwise.

#+BEGIN_SRC lisp
  (mlower-triangular-p m) => boolean
#+END_SRC
Return ~t~ if the matrix /M/ is in lower triangular form, and ~nil~ otherwise.

#+BEGIN_SRC lisp
  (mhessenberg-p m) => boolean
#+END_SRC
Return ~t~ if the matrix /M/ is in upper Hessenberg form, and ~nil~ otherwise.

*** Basic operations
#+BEGIN_SRC lisp
  (m+ m &rest matrices) => matrix
#+END_SRC
Return the addition of /matrices/ to the matrix /M/.

#+BEGIN_SRC lisp
  (m- m &rest matrices) => matrix
#+END_SRC
Return the subtraction of /matrices/ from the matrix /M/.

#+BEGIN_SRC lisp
  (m* m &rest matrices) => matrix
#+END_SRC
Return the product of the matrix /M/ and /matrices/.

#+BEGIN_SRC lisp
  (mc+ m x) => matrix
#+END_SRC
Return a matrix with the all elements of /M/ incremented by a constant /x/.

#+BEGIN_SRC lisp
  (mc- m x) => matrix
#+END_SRC
Return a matrix with the all elements of /M/ decremented by a constant /x/.

#+BEGIN_SRC lisp
  (mc* m x) => matrix
#+END_SRC
Return the mutiplication of the matrix /M/ by a constant /x/.

#+BEGIN_SRC lisp
  (mc/ m x) => matrix
#+END_SRC
Return the division of the matrix /M/ by a constant /x/.

#+BEGIN_SRC lisp
  (me* m &rest matrices) => matrix
#+END_SRC
Return the Hadamard (element-wise) product of the matrix /M/ and /matrices/.

#+BEGIN_SRC lisp
  (me/ m &rest matrices) => matrix
#+END_SRC
Return the element-wise division of the matrix /M/ by /matrices/.

#+BEGIN_SRC lisp
  (mk* m &rest matrices) => matrix
#+END_SRC
Return the Kronecker product of the matrix /M/ and /matrices/.

#+BEGIN_SRC lisp
  (m.* m1 m2) => number
#+END_SRC
Return the scalar product of the column matrices /M1/ and /M2/.

#+BEGIN_SRC lisp
  (mnorm m) => number
#+END_SRC
Return the Frobenius norm of the matrix /M/.

#+BEGIN_SRC lisp
  (mconj m) => matrix
#+END_SRC
Return the complex conjugate of the matrix /M/.

#+BEGIN_SRC lisp
  (mt m) => matrix
#+END_SRC
Return the transpose of the matrix /M/.

#+BEGIN_SRC lisp
  (mct m) => matrix
#+END_SRC
Return the conjugate transpose of the matrix /M/.

#+BEGIN_SRC lisp
  (mtr m) => number
#+END_SRC
Return the trace of the matrix /M/.

#+BEGIN_SRC lisp
  (minv m) => matrix, number
#+END_SRC
Return the inverse of the matrix /M/ and the determinant of /M/.
If /M/ is not invertible, ~nil~ and 0 will be returned.

#+BEGIN_SRC lisp
  (mpinv m) => matrix
#+END_SRC
Return the Moore-Penrose pseudoinverse of the matrix /M/.

#+BEGIN_SRC lisp
  (mdet m) => number
#+END_SRC
Return the determinant of the matrix /M/.

#+BEGIN_SRC lisp
  (mrk m) => number
#+END_SRC
Return the rank of the matrix /M/.

#+BEGIN_SRC lisp
  (mcom m) => matrix
#+END_SRC
Return the cofactor matrix of /M/.

#+BEGIN_SRC lisp
  (mexpt m n) => matrix
#+END_SRC
Return the matrix /M/ raised to the nth power.

#+BEGIN_SRC lisp
  (mexp m &optional iterations diagonalize-p) => matrix
#+END_SRC
Return the exponential of the matrix /M/. If /diagonalize-p/ is ~t~ (the
default), the function will try to diagonalize /M/ to compute its exponential.
If it fails, or if /diagonalize-p/ is ~nil~, the Taylor series will be used
instead.

#+BEGIN_SRC lisp
  (msqrt m &optional iterations diagonalize-p) => matrix
#+END_SRC
Return a square root of the matrix /M/. If /diagonalize-p/ is ~t~ (the
default), the function will try to diagonalize /M/ to compute its square root.
If it fails, or if /diagonalize-p/ is ~nil~, Heron's method will be used
instead. If no square root has been found, ~nil~ is returned.

#+BEGIN_SRC lisp
  (mfun m f &optional iterations) => matrix
#+END_SRC
Apply the matrix function /f/ to the matrix /M/. /f/ must take a number as
argument and return a number. If /M/ can't be diagonalized to apply /f/,
~nil~ is returned.

#+BEGIN_SRC lisp
  (mefun m f) => matrix
#+END_SRC
Return a matrix created by applying the function /f/ to each element of the
matrix /M/. /f/ must return a number, and take 3 numbers as arguments:
 - the element of /M/
 - the row index of the element
 - the column index of the element.

#+BEGIN_SRC lisp
  (mjacobian f x &optional dx) => matrix
#+END_SRC
Return the Jacobian matrix for the vector-valued function /f/ evaluated at /X/.
/f/ must take a column matrix as argument and return a column matrix. A step
/DX/ will be used to compute the finite differences. It can be a column matrix
with the same size as /X/ or a real number.

*** Decompositions
#+BEGIN_SRC lisp
  (mqr m &optional variant) => matrix, matrix
#+END_SRC
Return two matrices /Q/ and /R/, where /Q/ is orthonormal, /R/ is upper
triangular and /M = QR/.
If /variant/ is ~:householder~ or ~:givens~, a full QR decomposition is
returned.
If /variant/ is ~:householder-reduced~, ~:givens-reduced~ or ~:gram-schmidt~,
a reduced QR decomposition is returned.

#+BEGIN_SRC lisp
  (mlq m) => matrix, matrix
#+END_SRC
Return two matrices /L/ and /Q/, where /L/ is lower triangular, /Q/ is
orthonormal and /M = LQ/.

#+BEGIN_SRC lisp
  (mrq m) => matrix, matrix
#+END_SRC
Return two matrices /R/ and /Q/, where /R/ is upper triangular, /Q/ is
orthonormal, and /M = RQ/.

#+BEGIN_SRC lisp
  (mql m) => matrix, matrix
#+END_SRC
Return two matrices /Q/ and /L/, where /Q/ is orthonormal, /L/ is lower
triangular and /M = QL/.

#+BEGIN_SRC lisp
  (mlu m) => matrix, matrix
#+END_SRC
Return two matrices /L/ and /U/, where /L/ is lower triangular, /U/ is upper
triangular and /M = LU/.
If a LU decomposition is not possible, ~nil~ and ~nil~ are returned.

#+BEGIN_SRC lisp
  (mlup m) => matrix, matrix, matrix
#+END_SRC
Return 3 matrices /L/ and /U/ and /P/, where /L/ is lower triangular, /U/ is
upper triangular, /P/ is a permutation matrix and /PM = LU/.

#+BEGIN_SRC lisp
  (mcholesky m) => matrix
#+END_SRC
Return a lower triangular matrix /L/ such that /M = LL*/. /M/ must be
a hermitian positive-definite matrix. If a Cholesky decomposition is not
possible, ~nil~ is returned.

#+BEGIN_SRC lisp
  (mldl m) => matrix, matrix
#+END_SRC
Return 2 matrices /L/ and /D/, where /L/ is lower triangular with ones on the
diagonal, /D/ is diagonal, and /M = LDL*/. /M/ must be a hermitian matrix. If
an alternative Cholesky decomposition is not possible, ~nil~ and ~nil~ are
returned.

#+BEGIN_SRC lisp
  (mhessenberg m) => matrix, matrix
#+END_SRC
Return 2 matrices /H/ and /Q/, where /H/ is in upper Hessenberg form, /Q/ is
unitary, and /M = QHQ*/."

#+BEGIN_SRC lisp
  (mbalance m) => matrix, matrix
#+END_SRC
Return 2 matrices /A/ and /D/, such that /A/ is a balanced matrix, /D/ is
diagonal, and /M = DAD⁻¹/.

#+BEGIN_SRC lisp
  (mschur m &optional iterations variant) => matrix, matrix
#+END_SRC
Return 2 matrices /U/ and /Q/, where /U/ is upper triangular with the
eigenvalues on the diagonal, /Q/ is unitary, and /M = QUQ*/. If /variant/ is
~:real~ instead of ~:complex~, /U/ will be quasi-triangular if /M/ has complex
eigenvalues.

#+BEGIN_SRC lisp
  (mschur-reordering u indexes) => matrix, matrix
#+END_SRC
Given an upper triangular matrix /U/, return 2 matrices /V/ and /Q/ such that
/V/ is upper triangular, /Q/ is unitary and /U = QVQ*/ with the top left block
of /V/ containing the eigenvalues that are at some /indexes/ in the diagonal of
/U/.

#+BEGIN_SRC lisp
  (meigen m &optional iterations) => matrix, matrix
#+END_SRC
Return a diagonal matrix /D/ with the eigenvalues of the matrix /M/ on the
diagonal, and a matrix /P/ whose columns are eigenvectors. If /P/ is
invertible, /M = PDP⁻¹/.

#+BEGIN_SRC lisp
  (meigenvector u a &optional b iterations) => matrix
#+END_SRC
Return an eigenvector for the eigenvalue /u/. The returned column matrix /X/ is
such that /AX = uX/ if /B/ is not specified, or such that /AX = uBX/ if /B/ is
specified.

#+BEGIN_SRC lisp
  (mqz a b &optional iterations) => matrix, matrix, matrix, matrix
#+END_SRC
Return 4 matrices /R/, /S/, /Q/ and /Z/, where /R/ and /S/ are upper
triangular, /Q/ and /Z/ are unitary, /A = QRZ*/ and /B = QSZ*/.

#+BEGIN_SRC lisp
  (mqz-reordering r s indexes) => matrix, matrix, matrix, matrix
#+END_SRC
Given 2 matrices /R/ and /S/ in upper triangular form, return 4 matrices /U/,
/V/, /Q/ and /Z/ such that /U/ and /V/ are upper triangular, /Q/ and /Z/ are
unitary, /R = QUZ*/ and /S = QVZ*/ with the top left blocks of /U/ and /V/
containing the generalized eigenvalues that are at some /indexes/ in the
diagonals of /R/ and /S/.

#+BEGIN_SRC lisp
  (msvd m &optional iterations) => matrix, matrix, matrix
#+END_SRC
Return 3 matrices /U/, /S/ and /V/, where /U/ and /V/ are orthonormal,
/S/ is diagonal with the singular values on the diagonal, and /M = USV*/.

*** Solving equations
#+BEGIN_SRC lisp
  (msolve m y) => matrix, boolean
#+END_SRC
Return a column matrix /X/ such that /Y = MX/. The second returned value is ~t~
if /X/ is the unique solution, and ~nil~ otherwise. If no solution has been
found, ~nil~ and ~nil~ are returned.

#+BEGIN_SRC lisp
  (msolve-ode m x0 y u) => matrix
#+END_SRC
Return the column matrix /X(u)/, knowing that /dX(t)/dt = MX(t) + Y/
and /X(0) = X0/. If no solution has been found, ~nil~ is returned.

#+BEGIN_SRC lisp
  (msolve-axb y &rest matrices) => matrix, boolean
#+END_SRC
Given a matrix /Y/ and square /matrices/ /A/, /B/, /C/, /D/, /E/, /F/, ...,
return a matrix /X/ such that /AXB + CXD + EXF + ... = Y/. The second returned
value is ~t~ if /X/ is the unique solution, and ~nil~ otherwise. If no solution
has been found, ~nil~ and ~nil~ are returned.

#+BEGIN_SRC lisp
  (msolve-sylvester a b c) => matrix, boolean
#+END_SRC
Return a matrix /X/ such that /AX + XB = C/. The second returned value is ~t~
if /X/ is the unique solution, and ~nil~ otherwise. If no solution has been
found, ~nil~ and ~nil~ are returned.

#+BEGIN_SRC lisp
  (msolve-riccati a b c) => matrix
#+END_SRC
Return a matrix /X/ such that /AᵀX + XA - XBX + C = 0/.
If no solution has been found, ~nil~ is returned.

#+BEGIN_SRC lisp
  (msolve-riccati-discrete a b q r) => matrix
#+END_SRC
Return a matrix /X/ such that /AᵀXA - X - AᵀXB(R + BᵀXB)⁻¹BᵀXA + Q = 0/.
If no solution has been found, ~nil~ is returned.

#+BEGIN_SRC lisp
  (mls m y) => matrix
#+END_SRC
Return a column matrix /X/ minimizing the norm of /Y - MX/ (least squares).

#+BEGIN_SRC lisp
  (mls-axb y &rest matrices) => matrix
#+END_SRC
Given a matrix /Y/ and /matrices/ /A/, /B/, /C/, /D/, /E/, /F/, ...,
return a matrix /X/ minimizing the norm of /AXB + CXD + EXF + ... - Y/
(least squares).

* Examples
** Compute a number of the Fibonacci sequence
#+BEGIN_SRC lisp
  (asdf:load-system "simple-matrix")
  (use-package :simple-matrix)

  (defun fibonacci (n)
    (let* ((m (matrix 2 2 '(1 1 1 0)))
           (mn (mexpt m n)))
      (aref mn 0 1)))

  (fibonacci 39)

  => 63245986
#+END_SRC
We get: F(39) = 63245986

** System of equations with one exact solution
#+BEGIN_EXAMPLE
  2x + 3y = 10
  -x + 5y = 5
#+END_EXAMPLE
#+BEGIN_SRC lisp
  (let ((m (matrix 2 2 '(2 3 -1 5)))
        (y (matrix 2 1 '(10 5))))
    (mat->vec (msolve m y)))

  => #(35/13 20/13)
#+END_SRC
We get: x = 35/13, y = 20/13

** System of equations with an infinite number of solutions
One of the possible solutions will be returned.
#+BEGIN_EXAMPLE
  x + y + z = 6
  x - y + 3z = 4
  2x + 4z = 10
#+END_EXAMPLE
#+BEGIN_SRC lisp
  (let ((m (matrix 3 3 '(1 1 1 1 -1 3 2 0 4)))
        (y (matrix 3 1 '(6 4 10))))
    (mat->vec (mround (msolve m y))))

  => #(3.0d0 2.0d0 1.0d0)
#+END_SRC
We get: x = 3, y = 2, z = 1

** System of equations with no exact solution
We find an approximate solution using the method of least squares.
#+BEGIN_EXAMPLE
  2x + 3y = 10
  -x + 5y = 5
  3x - 3y = 7
#+END_EXAMPLE
#+BEGIN_SRC lisp
  (let ((m (matrix 3 2 '(2 3 -1 5 3 -3)))
        (y (matrix 3 1 '(10 5 7))))
    (mat->vec (mls m y)))

  => #(3.3828996282527886d0 1.42007434944238d0)
#+END_SRC
We get: x ≈ 3.383, y ≈ 1.42

[[./no-exact-solution.png]]

** Polynomial regression
Using the method of least squares, we find an approximation of the form:
#+BEGIN_EXAMPLE
y = c₀ + c₁x + c₂x² + c₃x³ + ... + cₙxⁿ
#+END_EXAMPLE
#+BEGIN_EXAMPLE
   X | Y
 ----|----
   2 |  1
   3 |  3
   6 | 17
   7 | 23
#+END_EXAMPLE
#+BEGIN_SRC lisp
  (defun polynomial-regression (degree x-values y-values)
    (let* ((n (length x-values))
           (y (vec->mat y-values))
           (m (matrix n (1+ degree))))
      ;; Fill the rows of the matrix with 1, x, x², ...
      (dotimes (i n)
        (dotimes (j (1+ degree))
          (setf (aref m i j) (expt (aref x-values i) j))))
      (mat->vec (mls m y))))

  (let ((x-values #(2 3 6 7))
        (y-values #(1 3 17 23)))
    (polynomial-regression 2 x-values y-values))

  => #(-1.117647058823536d0 -0.029411764705878803d0 0.49999999999999967d0)
#+END_SRC
We get the approximation:
#+BEGIN_EXAMPLE
y ≈ -1.118 - 0.029x + 0.5x²
#+END_EXAMPLE
[[./polynomial-regression.png]]

** Gaussian regression
Using the method of least squares, we find an approximation of the form:
#+BEGIN_EXAMPLE
y = A/(𝜎 * √(2𝜋)) * exp(-1/2 * ((x - 𝜇) / 𝜎)²)
#+END_EXAMPLE
Which can be transformed to:
#+BEGIN_EXAMPLE
ln(y) = c₀ + c₁x + c₂x²
#+END_EXAMPLE
With:
#+BEGIN_EXAMPLE
c₀ = ln(A / (𝜎 * √(2𝜋))) - 𝜇² / (2 * 𝜎²)
c₁ = 𝜇 / 𝜎²
c₂ = -1 / (2 * 𝜎²)

𝜎 = 1 / √(-2 * c₂)
𝜇 = c₁ * 𝜎²
A = 𝜎 * √(2𝜋) * exp(c₀ + 𝜇² / (2 * 𝜎²))
#+END_EXAMPLE
#+BEGIN_EXAMPLE
       X     |     Y
 ------------|------------
   5.647237  | 2.4841535
   3.8547695 | 1.0068833
   6.5579185 | 0.3458975
   5.8500857 | 1.8678124
   3.769868  | 0.82293093
   7.1886773 | 0.03394338
   6.6337585 | 0.28119978
   4.71034   | 3.758526
   3.8235922 | 1.0181205
   5.5816703 | 2.8982043
#+END_EXAMPLE
#+BEGIN_SRC lisp
  (defun gaussian-regression (x-values y-values)
    (let* ((n (length x-values))
           (y (vec->mat (map 'vector #'log y-values)))
           (m (matrix n 3)))
      (dotimes (i n)
        (dotimes (j 3)
          (setf (aref m i j) (expt (aref x-values i) j))))
      (let* ((x (mls m y))
             (c0 (aref x 0 0))
             (c1 (aref x 1 0))
             (c2 (aref x 2 0))
             (sigma (/ (sqrt (* -2 c2))))
             (mu (* c1 sigma sigma))
             (a (* sigma
                   (sqrt (* 2 pi))
                   (exp (+ c0 (/ (* mu mu) (* 2 sigma sigma)))))))
        (vector a mu sigma))))

  (let ((x-values #(5.647237 3.8547695 6.5579185 5.8500857 3.769868
                    7.1886773 6.6337585 4.71034 3.8235922 5.5816703))
        (y-values #(2.4841535 1.0068833 0.3458975 1.8678124 0.82293093
                    0.03394338 0.28119978 3.758526 1.0181205 2.8982043)))
    (gaussian-regression x-values y-values))

  => #(6.934017483512659d0 5.005451177669433d0 0.7077270060667742d0)
#+END_SRC
We get: A ≈ 6.934, 𝜇 ≈ 5.005 , 𝜎 ≈ 0.708; which gives the approximation:
#+BEGIN_EXAMPLE
y ≈ 3.909 * exp(-0.998 * (x - 5.005)²)
#+END_EXAMPLE
[[./gaussian-regression.png]]

** Markov chain
We find the stationary distribution of a system with 3 states, whose transition
matrix is:
#+BEGIN_EXAMPLE
⎛ 90/100  5/100   5/100 ⎞
⎜ 70/100    0    30/100 ⎟
⎝ 80/100    0    20/100 ⎠
#+END_EXAMPLE

[[./markov-chain.png]]

#+BEGIN_SRC lisp
  (defun markov-stationary (transition-matrix)
    (let* ((n (array-dimension transition-matrix 0))
           (m (m- transition-matrix (mi n)))
           (y (matrix 1 n)))
      ;; Replace a column of M and Y to indicate that x₀ + x₁ + ... + xₙ = 1
      (setf (aref y 0 0) 1)
      (dotimes (i n)
        (setf (aref m i 0) 1))
      (mat->vec (msolve (mt m) (mt y)))))

  (let ((m (matrix 3 3 '(90/100 5/100 5/100
                         70/100 0 30/100
                         80/100 0 20/100))))
    (markov-stationary m))

  => #(160/181 8/181 13/181)
#+END_SRC
We get:
 - probability to be in the first state: 160/181 ≈ 0.884
 - probability to be in the second state: 8/181 ≈ 0.044
 - probability to be in the third state: 13/181 ≈ 0.072

** System of differential equations
#+BEGIN_EXAMPLE
  df(x)/dx = f(x) - 2g(x) + 3
  dg(x)/dx = 3f(x) - 4g(x) - 2

  f(0) = 11
  g(0) = 19/2
#+END_EXAMPLE
#+BEGIN_SRC lisp
  (let ((m (matrix 2 2 '(1 -2 3 -4)))
        (v (matrix 2 1 '(11 19/2)))
        (w (matrix 2 1 '(3 -2)))
        (x 3))
    (mat->vec (mround (msolve-ode m v w x))))

  => #(8.054745d0 5.557223d0)
#+END_SRC
We get: f(3) ≈ 8.055, g(3) ≈ 5.557

[[./differential-equations.png]]

** Linear-quadratic regulator
We find the optimal control gain matrix for a linear-quadratic regulator for
a linear system defined by:
#+BEGIN_EXAMPLE
dx(t)/dx = Ax(t) + Bu(t)
y(t) = x(t)
u(t) = -Kx(t)

A = ⎛  0   1 ⎞   B = ⎛ 0 ⎞   Q = ⎛ 1 0 ⎞   R = ( 1 )
    ⎝ -2  -3 ⎠       ⎝ 1 ⎠       ⎝ 0 0 ⎠
#+END_EXAMPLE

[[./linear-system.png]]

#+BEGIN_SRC lisp
  (defun lqr (a b q r)
    (let* ((tb (mt b))
           (invr (minv r))
           (p (msolve-riccati a (m* b invr tb) q)))
      (m* invr tb p)))

  (let ((a (matrix 2 2 '(0 1 -2 -3)))
        (b (matrix 2 1 '(0 1)))
        (q (matrix 2 2 '(1 0 0 0)))
        (r (matrix 1 1 '(1))))
    (mround (lqr a b q r)))

  => #2A((0.236068d0 0.077684d0))
#+END_SRC
We get: K ≈ ( 0.236 0.078 )

* Tests
The tests require the [[https://common-lisp.net/project/fiveam/][fiveam]] package.
They can be run with:

#+BEGIN_SRC lisp
(asdf:test-system "simple-matrix")
#+END_SRC