Linear algebra package

Matrix decomposition

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:

Here are examples of its use:
  matrix_t T = Decomp<LU>(M); 
  matrix_t R = Decomp<QR>(M);
See also:

Solving a linear system

Several methods for solving linear systems

\[ A X = B \]

where $A$ and $B$ are $n \times m$ and $n \times p$ matrices. The result is a $m \times p$ matrix. On can use the following methods

Here are examplesof its use.
  matrix_t B = ...;
  matrix_t X = solve<LU>(A,B);
  matrix_t Y = solve<LS>(A,B);
See also:


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

See also:


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:


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

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

See also:

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:

Fast Fourier Transform

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

See also: