The general interface for matrix decomposition is given by the function:

template<class Mthd, class Mat> Mat decomp(const Mat & m);

The first parameter specifies the type of decomposition. It can be:

LU for Lower-Upper triangular decomposition.

Bareiss for Lower-Upper triangular decomposition based on Bareiss method. It applies the following pivoting scheme:
$a^{(k)}_{i, j} = \frac{a^{(k - 1)}_{i, j} a^{(k - 1)}_k - a^{(k - 1)}_{i, k} a^{(k - 1)}_{k, j}}{a^{(k - 2)}_{k - 1, k - 1}},$
where $A^{(k)} = (a_{i, j}^{(k)})$ is the transformed matrix at step and $A^{(0)} \assign A$ is the initial matrix. At the end of this process, after some permutation of the rows, the matrix is upper triangular. If A is a square matrix, the last entry $A^{(n)}_{n, n}$ on the diagonal at the end of the process is the determinant of the matrix . Moreover, if has its coefficients in a ring, it is the same for the matrices $A^{(k)}$ , since the previous division is exact

QR for Orthogonal-Upper triangular decomposition.

Here are examples of its use:

  matrix_t T = Decomp<LU>(M); 
  matrix_t R = Decomp<QR>(M);

See also:: synaps/linalg/decomp.h

Determinant

Different methods for computing the determinant of a matrix are discribed here: Gauss, Bareiss or Berkowitz methods.

LU for Gauss method uses classical pivoting techniques (ie. LU decomposition) to compute the determinant.

Bareiss for Bareiss method .

Berkowitz method is based on the following recursion formula, based on Cayley-Hamilton identity:
$d_0 : = \tmop{Trace} (A) ; A_0 : = A ; A_{k + 1} : = A_k A - \frac{1}{k + 1} \tmop{Trace} (A_k) A.$
A mixed strategy, which expands the determinant with respect to the first row and Bareiss method is applied to the subdeterminants.

See also:: synaps/linalg/Det.h

SVD

Theorem: [GLVL96] For a real matrix $A$

of size $m \times n$ , there exists two orthogonal matrices

$U = [u_1 \cdots, u_m] \in \mathbb{R}^{m \times m}, \tmop{and} V = [v_1, \cdots, \quad v_n] \in \mathbb{R}^{n \times n}$

such that

$U^t AV = \tmop{diag} (\sigma_1, \cdots, \sigma_p) \in \mathbb{R}^{m \times n} \quad p = \min \{m, n\} \tmtextrm{,}$

$\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_p \geq 0 \hspace{2em} \hspace{2em} \hspace{2em} \hspace{2em} \hspace{2em} \hspace{2em}$

$\sigma_i$ are called the $i$ th singular value of A, $u_i$ and $v_i$ are respectively the $i$ th left and right singular vectors.

Definition: For a matrix $A$ , the numerical rank of A is defined by:
$\tmop{rank} (A, \epsilon) = \min \{\tmop{rank} (B) : \left| \frac{A}{\|A\|} - B \right| \leq \epsilon\}.$

Theorem: [GLVL96] If $A \in \mathbb{R}^{m \times n}$ is of rank $r$ and $k < r$ then
$\sigma_{k + 1} = \min_{\tmop{rank} (B) = k} \|A - B\|.$

The singular value decomposition of a matrix $M$ is the decomposition of the form

$M = U \Sigma V^t$

where $U$ and $V$ are orthognal matrices $U U^t = \tmop{Id}$ , $V V^t = \tmop{Id}$ and $\Sigma = \tmop{diag} (\sigma_1, \ldots, \sigma_n)$ is a diagonal matrix such that $\sigma_1 \geqslant \cdots \geqslant \sigma_n \geqslant 0$ . The values $\sigma_i$ are called the singular values of $M$ .

  vector_t s =  Svd(m);

See also:: synaps/linalg/svd.h

Rank

Theorem [*] shows that the smallest singular value is the distance between A and all the matrices of rank $< p = \min \{m, n\}$ . Thus if $r_{\epsilon} = \tmop{rank} (A, \epsilon)$ then

$\sigma_1 \geq \cdots \geq \sigma_{r_{\epsilon}} > \epsilon \sigma_1 \geq \sigma_{r_{\epsilon + 1}} \geq \cdots \geq \sigma_p, \quad p = \min \{m, n\}.$

As a result, in order to determine the rank of a matrix $A$ , we find the singular values $(\sigma_i)_i$ so that

$\tmop{rank} (A, \epsilon) = \max \{i ; \frac{\sigma_i}{\sigma_1} > \epsilon\}.$

The function

  Rank(A);

computes the rank of a matrix, as described above with a default value for $\varepsilon$ .

See also:: synaps/linalg/rank.h

Eigenvalues, eigenvectors

The computation of the eigenvalues and eigenvectors of a matrix are performed by the functions eigen and eigen_vect. They also provide generalised eigenvalues and eigenvectors, for pairs of matrices (a,B). If the internal representation of type lapack::rep2d<double>, the routine dgeeg or dgegv are called.

  VectDse<std::complex<double> > v = eigen(A);
  VectDse<double>                w = eigen(A, Real());
  VectDse<double>                w = eigen(A,B);

See also:: synaps/linalg/eigen.h

Fast Fourier Transform

A set of functions are available, for performing FFT on generic vectors.

See also:: synaps/linalg/FFT.h

SYNAPS DOCUMENTATION

Linear algebra package

Matrix decomposition

Solving a linear system

Determinant

SVD

Rank

Eigenvalues, eigenvectors

Fast Fourier Transform