// Copyright (C) 2009 Davis E. King (davis@dlib.net)
// License: Boost Software License See LICENSE.txt for the full license.
#undef DLIB_MATRIx_LA_FUNCTS_ABSTRACT_
#ifdef DLIB_MATRIx_LA_FUNCTS_ABSTRACT_
#include "matrix_abstract.h"
#include <complex>
namespace dlib
{
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// Global linear algebra functions
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
const matrix_exp::matrix_type inv (
const matrix_exp& m
);
/*!
requires
- m is a square matrix
ensures
- returns the inverse of m
(Note that if m is singular or so close to being singular that there
is a lot of numerical error then the returned matrix will be bogus.
You can check by seeing if m*inv(m) is an identity matrix)
!*/
// ----------------------------------------------------------------------------------------
const matrix pinv (
const matrix_exp& m,
double tol = 0
);
/*!
requires
- tol >= 0
ensures
- returns the Moore-Penrose pseudoinverse of m.
- The returned matrix has m.nc() rows and m.nr() columns.
- if (tol == 0) then
- singular values less than max(m.nr(),m.nc()) times the machine epsilon
times the largest singular value are ignored.
- else
- singular values less than tol*max(singular value in m) are ignored.
!*/
// ----------------------------------------------------------------------------------------
void svd (
const matrix_exp& m,
matrix<matrix_exp::type>& u,
matrix<matrix_exp::type>& w,
matrix<matrix_exp::type>& v
);
/*!
ensures
- computes the singular value decomposition of m
- m == #u*#w*trans(#v)
- trans(#u)*#u == identity matrix
- trans(#v)*#v == identity matrix
- diag(#w) == the singular values of the matrix m in no
particular order. All non-diagonal elements of #w are
set to 0.
- #u.nr() == m.nr()
- #u.nc() == m.nc()
- #w.nr() == m.nc()
- #w.nc() == m.nc()
- #v.nr() == m.nc()
- #v.nc() == m.nc()
- if DLIB_USE_LAPACK is #defined then the xGESVD routine
from LAPACK is used to compute the SVD.
!*/
// ----------------------------------------------------------------------------------------
long svd2 (
bool withu,
bool withv,
const matrix_exp& m,
matrix<matrix_exp::type>& u,
matrix<matrix_exp::type>& w,
matrix<matrix_exp::type>& v
);
/*!
requires
- m.nr() >= m.nc()
ensures
- computes the singular value decomposition of matrix m
- m == subm(#u,get_rect(m))*diagm(#w)*trans(#v)
- trans(#u)*#u == identity matrix
- trans(#v)*#v == identity matrix
- #w == the singular values of the matrix m in no
particular order.
- #u.nr() == m.nr()
- #u.nc() == m.nr()
- #w.nr() == m.nc()
- #w.nc() == 1
- #v.nr() == m.nc()
- #v.nc() == m.nc()
- if (widthu == false) then
- ignore the above regarding #u, it isn't computed and its
output state is undefined.
- if (widthv == false) then
- ignore the above regarding #v, it isn't computed and its
output state is undefined.
- returns an error code of 0, if no errors and 'k' if we fail to
converge at the 'kth' singular value.
- if (DLIB_USE_LAPACK is #defined) then
- if (withu == withv) then
- the xGESDD routine from LAPACK is used to compute the SVD.
- else
- the xGESVD routine from LAPACK is used to compute the SVD.
!*/
// ----------------------------------------------------------------------------------------
void svd3 (
const matrix_exp& m,
matrix<matrix_exp::type>& u,
matrix<matrix_exp::type>& w,
matrix<matrix_exp::type>& v
);
/*!
ensures
- computes the singular value decomposition of m
- m == #u*diagm(#w)*trans(#v)
- trans(#u)*#u == identity matrix
- trans(#v)*#v == identity matrix
- #w == the singular values of the matrix m in no
particular order.
- #u.nr() == m.nr()
- #u.nc() == m.nc()
- #w.nr() == m.nc()
- #w.nc() == 1
- #v.nr() == m.nc()
- #v.nc() == m.nc()
- if DLIB_USE_LAPACK is #defined then the xGESVD routine
from LAPACK is used to compute the SVD.
!*/
// ----------------------------------------------------------------------------------------
template <
typename T
>
void svd_fast (
const matrix<T>& A,
matrix<T>& u,
matrix<T>& w,
matrix<T>& v,
unsigned long l,
unsigned long q = 1
);
/*!
requires
- l > 0
- A.size() > 0
(i.e. A can't be an empty matrix)
ensures
- computes the singular value decomposition of A.
- Lets define some constants we use to document the behavior of svd_fast():
- Let m = A.nr()
- Let n = A.nc()
- Let k = min(l, min(m,n))
- Therefore, A represents an m by n matrix and svd_fast() is designed
to find a rank-k representation of it.
- if (the rank of A is <= k) then
- A == #u*diagm(#w)*trans(#v)
- else
- A is approximated by #u*diagm(#w)*trans(#v)
(i.e. In this case A can't be represented with a rank-k matrix, so the
matrix you get by trying to reconstruct A from the output of the SVD is
not exactly the same.)
- trans(#u)*#u == identity matrix
- trans(#v)*#v == identity matrix
- #w == the top k singular values of the matrix A (in no particular order).
- #u.nr() == m
- #u.nc() == k
- #w.nr() == k
- #w.nc() == 1
- #v.nr() == n
- #v.nc() == k
- This function implements the randomized subspace iteration defined in the
algorithm 4.4 and 5.1 boxes of the paper:
Finding Structure with Randomness: Probabilistic Algorithms for
Constructing Approximate Matrix Decompositions by Halko et al.
Therefore, it is very fast and suitable for use with very large matrices.
Moreover, q is the number of subspace iterations performed. Larger
values of q might increase the accuracy of the solution but the default
value should be good for many problems.
!*/
// ----------------------------------------------------------------------------------------
template <
typename sparse_vector_type,
typename T
>
void svd_fast (
const std::vector<sparse_vector_type>& A,
matrix<T>& u,
matrix<T>& w,
matrix<T>& v,
unsigned long l,
unsigned long q = 1
);
/*!
requires
- A contains a set of sparse vectors. See dlib/svm/sparse_vector_abstract.h
for a definition of what constitutes a sparse vector.
- l > 0
- max_index_plus_one(A) > 0
(i.e. A can't be an empty matrix)
ensures
- computes the singular value decomposition of A. In this case, we interpret A
as a matrix of A.size() rows, where each row is defined by a sparse vector.
- Lets define some constants we use to document the behavior of svd_fast():
- Let m = A.size()
- Let n = max_index_plus_one(A)
- Let k = min(l, min(m,n))
- Therefore, A represents an m by n matrix and svd_fast() is designed
to find a rank-k representation of it.
- if (the rank of A is <= k) then
- A == #u*diagm(#w)*trans(#v)
- else
- A is approximated by #u*diagm(#w)*trans(#v)
(i.e. In this case A can't be represented with a rank-k matrix, so the
matrix you get by trying to reconstruct A from the output of the SVD is
not exactly the same.)
- trans(#u)*#u == identity matrix
- trans(#v)*#v == identity matrix
- #w == the top k singular values of the matrix A (in no particular order).
- #u.nr() == m
- #u.nc() == k
- #w.nr() == k
- #w.nc() == 1
- #v.nr() == n
- #v.nc() == k
- This function implements the randomized subspace iteration defined in the
algorithm 4.4 and 5.1 boxes of the paper:
Finding Structure with Randomness: Probabilistic Algorithms for
Constructing Approximate Matrix Decompositions by Halko et al.
Therefore, it is very fast and suitable for use with very large matrices.
Moreover, q is the number of subspace iterations performed. Larger
values of q might increase the accuracy of the solution but the default
value should be good for many problems.
!*/
// ----------------------------------------------------------------------------------------
template <
typename T,
long NR,
long NC,
typename MM,
typename L
>
void orthogonalize (
matrix<T,NR,NC,MM,L>& m
);
/*!
requires
- m.nr() >= m.nc()
- m.size() > 0
ensures
- #m == an orthogonal matrix with the same dimensions as m. In particular,
the columns of #m have the same span as the columns of m.
- trans(#m)*#m == identity matrix
- This function is just shorthand for computing the QR decomposition of m
and then storing the Q factor into #m.
!*/
// ----------------------------------------------------------------------------------------
const matrix real_eigenvalues (
const matrix_exp& m
);
/*!
requires
- m.nr() == m.nc()
- matrix_exp::type == float or double
ensures
- returns a matrix E such that:
- E.nr() == m.nr()
- E.nc() == 1
- E contains the real part of all eigenvalues of the matrix m.
(note that the eigenvalues are not sorted)
!*/
// ----------------------------------------------------------------------------------------
const matrix_exp::type det (
const matrix_exp& m
);
/*!
requires
- m is a square matrix
ensures
- returns the determinant of m
!*/
// ----------------------------------------------------------------------------------------
const matrix_exp::type trace (
const matrix_exp& m
);
/*!
requires
- m is a square matrix
ensures
- returns the trace of m
(i.e. returns sum(diag(m)))
!*/
// ----------------------------------------------------------------------------------------
const matrix_exp::matrix_type chol (
const matrix_exp& A
);
/*!
requires
- A is a square matrix
ensures
- if (A has a Cholesky Decomposition) then
- returns the decomposition of A. That is, returns a matrix L
such that L*trans(L) == A. L will also be lower triangular.
- else
- returns a matrix with the same dimensions as A but it
will have a bogus value. I.e. it won't be a decomposition.
In this case the algorithm returns a partial decomposition.
- You can tell when chol fails by looking at the lower right
element of the returned matrix. If it is 0 then it means
A does not have a cholesky decomposition.
- If DLIB_USE_LAPACK is defined then the LAPACK routine xPOTRF
is used to compute the cholesky decomposition.
!*/
// ----------------------------------------------------------------------------------------
const matrix_exp::matrix_type inv_lower_triangular (
const matrix_exp& A
);
/*!
requires
- A is a square matrix
ensures
- if (A is lower triangular) then
- returns the inverse of A.
- else
- returns a matrix with the same dimensions as A but it
will have a bogus value. I.e. it won't be an inverse.
!*/
// ----------------------------------------------------------------------------------------
const matrix_exp::matrix_type inv_upper_triangular (
const matrix_exp& A
);
/*!
requires
- A is a square matrix
ensures
- if (A is upper triangular) then
- returns the inverse of A.
- else
- returns a matrix with the same dimensions as A but it
will have a bogus value. I.e. it won't be an inverse.
!*/
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
// Matrix decomposition classes
// ----------------------------------------------------------------------------------------
// ----------------------------------------------------------------------------------------
template <
typename matrix_exp_type
>
class lu_decomposition
{
/*!
REQUIREMENTS ON matrix_exp_type
must be some kind of matrix expression as defined in the
dlib/matrix/matrix_abstract.h file. (e.g. a dlib::matrix object)
The matrix type must also contain float or double values.
WHAT THIS OBJECT REPRESENTS
This object represents something that can compute an LU
decomposition of a real valued matrix. That is, for any
matrix A it computes matrices L, U, and a pivot vector P such
that rowm(A,P) == L*U.
The LU decomposition with pivoting always exists, even if the matrix is
singular, so the constructor will never fail. The primary use of the
LU decomposition is in the solution of square systems of simultaneous
linear equations. This will fail if is_singular() returns true (or
if A is very nearly singular).
If DLIB_USE_LAPACK is defined then the LAPACK routine xGETRF
is used to compute the LU decomposition.
!*/
public:
const static long NR = matrix_exp_type::NR;
const static long NC = matrix_exp_type::NC;
typedef typename matrix_exp_type::type type;
typedef typename matrix_exp_type::mem_manager_type mem_manager_type;
typedef typename matrix_exp_type::layout_type layout_type;
typedef matrix<type,0,0,mem_manager_type,layout_type> matrix_type;
typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type;
typedef matrix<long,NR,1,mem_manager_type,layout_type> pivot_column_vector_type;
template <typename EXP>
lu_decomposition (
const matrix_exp<EXP> &A
);
/*!
requires
- EXP::type == lu_decomposition::type
- A.size() > 0
ensures
- #nr() == A.nr()
- #nc() == A.nc()
- #is_square() == (A.nr() == A.nc())
- computes the LU factorization of the given A matrix.
!*/
bool is_square (
) const;
/*!
ensures
- if (the input A matrix was a square matrix) then
- returns true
- else
- returns false
!*/
bool is_singular (
) const;
/*!
requires
- is_square() == true
ensures
- if (the input A matrix is singular) then
- returns true
- else
- returns false
!*/
long nr(
) const;
/*!
ensures
- returns the number of rows in the input matrix
!*/
long nc(
) const;
/*!
ensures
- returns the number of columns in the input matrix
!*/
const matrix_type get_l (
) const;
/*!
ensures
- returns the lower triangular L factor of the LU factorization.
- L.nr() == nr()
- L.nc() == min(nr(),nc())
!*/
const matrix_type get_u (
) const;
/*!
ensures
- returns the upper triangular U factor of the LU factorization.
- U.nr() == min(nr(),nc())
- U.nc() == nc()
!*/
const pivot_column_vector_type& get_pivot (
) const;
/*!
ensures
- returns the pivot permutation vector. That is,
if A is the input matrix then this function
returns a vector P such that:
- rowm(A,P) == get_l()*get_u()
- P.nr() == A.nr()
!*/
type det (
) const;
/*!
requires
- is_square() == true
ensures
- computes and returns the determinant of the input
matrix using LU factors.
!*/
template <typename EXP>
const matrix_type solve (
const matrix_exp<EXP> &B
) const;
/*!
requires
- EXP::type == lu_decomposition::type
- is_square() == true
- B.nr() == nr()
ensures
- Let A denote the input matrix to this class's constructor.
Then this function solves A*X == B for X and returns X.
- Note that if A is singular (or very close to singular) then
the X returned by this function won't fit A*X == B very well (if at all).
!*/
};
// ----------------------------------------------------------------------------------------
template <
typename matrix_exp_type
>
class cholesky_decomposition
{
/*!
REQUIREMENTS ON matrix_exp_type
must be some kind of matrix expression as defined in the
dlib/matrix/matrix_abstract.h file. (e.g. a dlib::matrix object)
The matrix type must also contain float or double values.
WHAT THIS OBJECT REPRESENTS
This object represents something that can compute a cholesky
decomposition of a real valued matrix. That is, for any
symmetric, positive definite matrix A, it computes a lower
triangular matrix L such that A == L*trans(L).
If the matrix is not symmetric or positive definite, the function
computes only a partial decomposition. This can be tested with
the is_spd() flag.
If DLIB_USE_LAPACK is defined then the LAPACK routine xPOTRF
is used to compute the cholesky decomposition.
!*/
public:
const static long NR = matrix_exp_type::NR;
const static long NC = matrix_exp_type::NC;
typedef typename matrix_exp_type::type type;
typedef typename matrix_exp_type::mem_manager_type mem_manager_type;
typedef typename matrix_exp_type::layout_type layout_type;
typedef typename matrix_exp_type::matrix_type matrix_type;
typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type;
template <typename EXP>
cholesky_decomposition(
const matrix_exp<EXP>& A
);
/*!
requires
- EXP::type == cholesky_decomposition::type
- A.size() > 0
- A.nr() == A.nc()
(i.e. A must be a square matrix)
ensures
- if (A is symmetric positive-definite) then
- #is_spd() == true
- Constructs a lower triangular matrix L, such that L*trans(L) == A.
and #get_l() == L
- else
- #is_spd() == false
!*/
bool is_spd(
) const;
/*!
ensures
- if (the input matrix was symmetric positive-definite) then
- returns true
- else
- returns false
!*/
const matrix_type& get_l(
) const;
/*!
ensures
- returns the lower triangular factor, L, such that L*trans(L) == A
(where A is the input matrix to this class's constructor)
- Note that if A is not symmetric positive definite or positive semi-definite
then the equation L*trans(L) == A won't hold.
!*/
template <typename EXP>
const matrix solve (
const matrix_exp<EXP>& B
) const;
/*!
requires
- EXP::type == cholesky_decomposition::type
- B.nr() == get_l().nr()
(i.e. the number of rows in B must match the number of rows in the
input matrix A)
ensures
- Let A denote the input matrix to this class's constructor. Then
this function solves A*X = B for X and returns X.
- Note that if is_spd() == false or A was really close to being
non-SPD then the solver will fail to find an accurate solution.
!*/
};
// ----------------------------------------------------------------------------------------
template <
typename matrix_exp_type
>
class qr_decomposition
{
/*!
REQUIREMENTS ON matrix_exp_type
must be some kind of matrix expression as defined in the
dlib/matrix/matrix_abstract.h file. (e.g. a dlib::matrix object)
The matrix type must also contain float or double values.
WHAT THIS OBJECT REPRESENTS
This object represents something that can compute a classical
QR decomposition of an m-by-n real valued matrix A with m >= n.
The QR decomposition is an m-by-n orthogonal matrix Q and an
n-by-n upper triangular matrix R so that A == Q*R. The QR decomposition
always exists, even if the matrix does not have full rank, so the
constructor will never fail. The primary use of the QR decomposition
is in the least squares solution of non-square systems of simultaneous
linear equations. This will fail if is_full_rank() returns false or
A is very nearly not full rank.
The Q and R factors can be retrieved via the get_q() and get_r()
methods. Furthermore, a solve() method is provided to find the
least squares solution of Ax=b using the QR factors.
If DLIB_USE_LAPACK is #defined then the xGEQRF routine
from LAPACK is used to compute the QR decomposition.
!*/
public:
const static long NR = matrix_exp_type::NR;
const static long NC = matrix_exp_type::NC;
typedef typename matrix_exp_type::type type;
typedef typename matrix_exp_type::mem_manager_type mem_manager_type;
typedef typename matrix_exp_type::layout_type layout_type;
typedef matrix<type,0,0,mem_manager_type,layout_type> matrix_type;
template <typename EXP>
qr_decomposition(
const matrix_exp<EXP>& A
);
/*!
requires
- EXP::type == qr_decomposition::type
- A.nr() >= A.nc()
- A.size() > 0
ensures
- #nr() == A.nr()
- #nc() == A.nc()
- computes the QR decomposition of the given A matrix.
!*/
bool is_full_rank(
) const;
/*!
ensures
- if (the input A matrix had full rank) then
- returns true
- else
- returns false
!*/
long nr(
) const;
/*!
ensures
- returns the number of rows in the input matrix
!*/
long nc(
) const;
/*!
ensures
- returns the number of columns in the input matrix
!*/
const matrix_type get_r (
) const;
/*!
ensures
- returns a matrix R such that:
- R is the upper triangular factor, R, of the QR factorization
- get_q()*R == input matrix A
- R.nr() == nc()
- R.nc() == nc()
!*/
const matrix_type get_q (
) const;
/*!
ensures
- returns a matrix Q such that:
- Q is the economy-sized orthogonal factor Q from the QR
factorization.
- trans(Q)*Q == identity matrix
- Q*get_r() == input matrix A
- Q.nr() == nr()
- Q.nc() == nc()
!*/
void get_q (
matrix_type& Q
) const;
/*!
ensures
- #Q == get_q()
- This function exists to allow a user to get the Q matrix without the
overhead of returning a matrix by value.
!*/
template <typename EXP>
const matrix_type solve (
const matrix_exp<EXP>& B
) const;
/*!
requires
- EXP::type == qr_decomposition::type
- B.nr() == nr()
ensures
- Let A denote the input matrix to this class's constructor.
Then this function finds the least squares solution to the equation A*X = B
and returns X. X has the following properties:
- X is the matrix that minimizes the two norm of A*X-B. That is, it
minimizes sum(squared(A*X - B)).
- X.nr() == nc()
- X.nc() == B.nc()
- Note that this function will fail to output a good solution if is_full_rank() == false
or the A matrix is close to not being full rank.
!*/
};
// ----------------------------------------------------------------------------------------
template <
typename matrix_exp_type
>
class eigenvalue_decomposition
{
/*!
REQUIREMENTS ON matrix_exp_type
must be some kind of matrix expression as defined in the
dlib/matrix/matrix_abstract.h file. (e.g. a dlib::matrix object)
The matrix type must also contain float or double values.
WHAT THIS OBJECT REPRESENTS
This object represents something that can compute an eigenvalue
decomposition of a real valued matrix. So it gives
you the set of eigenvalues and eigenvectors for a matrix.
Let A denote the input matrix to this object's constructor. Then
what this object does is it finds two matrices, D and V, such that
- A*V == V*D
Where V is a square matrix that contains all the eigenvectors
of the A matrix (each column of V is an eigenvector) and
D is a diagonal matrix containing the eigenvalues of A.
It is important to note that if A is symmetric or non-symmetric you
get somewhat different results. If A is a symmetric matrix (i.e. A == trans(A))
then:
- All the eigenvalues and eigenvectors of A are real numbers.
- Because of this there isn't really any point in using the
part of this class's interface that returns complex matrices.
All you need are the get_real_eigenvalues() and
get_pseudo_v() functions.
- V*trans(V) should be equal to the identity matrix. That is, all the
eigenvectors in V should be orthonormal.
- So A == V*D*trans(V)
- If DLIB_USE_LAPACK is #defined then this object uses the xSYEVR LAPACK
routine.
On the other hand, if A is not symmetric then:
- Some of the eigenvalues and eigenvectors might be complex numbers.
- An eigenvalue is complex if and only if its corresponding eigenvector
is complex. So you can check for this case by just checking
get_imag_eigenvalues() to see if any values are non-zero. You don't
have to check the V matrix as well.
- V*trans(V) won't be equal to the identity matrix but it is usually
invertible. So A == V*D*inv(V) is usually a valid statement but
A == V*D*trans(V) won't be.
- If DLIB_USE_LAPACK is #defined then this object uses the xGEEV LAPACK
routine.
!*/
public:
const static long NR = matrix_exp_type::NR;
const static long NC = matrix_exp_type::NC;
typedef typename matrix_exp_type::type type;
typedef typename matrix_exp_type::mem_manager_type mem_manager_type;
typedef typename matrix_exp_type::layout_type layout_type;
typedef typename matrix_exp_type::matrix_type matrix_type;
typedef matrix<type,NR,1,mem_manager_type,layout_type> column_vector_type;
typedef matrix<std::complex<type>,0,0,mem_manager_type,layout_type> complex_matrix_type;
typedef matrix<std::complex<type>,NR,1,mem_manager_type,layout_type> complex_column_vector_type;
template <typename EXP>
eigenvalue_decomposition(
const matrix_exp<EXP>& A
);
/*!
requires
- A.nr() == A.nc()
- A.size() > 0
- EXP::type == eigenvalue_decomposition::type
ensures
- #dim() == A.nr()
- computes the eigenvalue decomposition of A.
- #get_eigenvalues() == the eigenvalues of A
- #get_v() == all the eigenvectors of A
!*/
template <typename EXP>
eigenvalue_decomposition(
const matrix_op<op_make_symmetric<EXP> >& A
);
/*!
requires
- A.nr() == A.nc()
- A.size() > 0
- EXP::type == eigenvalue_decomposition::type
ensures
- #dim() == A.nr()
- computes the eigenvalue decomposition of the symmetric matrix A. Does so
using a method optimized for symmetric matrices.
- #get_eigenvalues() == the eigenvalues of A
- #get_v() == all the eigenvectors of A
- moreover, since A is symmetric there won't be any imaginary eigenvalues. So
we will have:
- #get_imag_eigenvalues() == 0
- #get_real_eigenvalues() == the eigenvalues of A
- #get_pseudo_v() == all the eigenvectors of A
- diagm(#get_real_eigenvalues()) == #get_pseudo_d()
Note that the symmetric matrix operator is created by the
dlib::make_symmetric() function. This function simply reflects
the lower triangular part of a square matrix into the upper triangular
part to create a symmetric matrix. It can also be used to denote that a
matrix is already symmetric using the C++ type system.
!*/
long dim (
) const;
/*!
ensures
- dim() == the number of rows/columns in the input matrix A
!*/
const complex_column_vector_type get_eigenvalues (
) const;
/*!
ensures
- returns diag(get_d()). That is, returns a
vector that contains the eigenvalues of the input
matrix.
- the returned vector has dim() rows
- the eigenvalues are not sorted in any particular way
!*/
const column_vector_type& get_real_eigenvalues (
) const;
/*!
ensures
- returns the real parts of the eigenvalues. That is,
returns real(get_eigenvalues())
- the returned vector has dim() rows
- the eigenvalues are not sorted in any particular way
!*/
const column_vector_type& get_imag_eigenvalues (
) const;
/*!
ensures
- returns the imaginary parts of the eigenvalues. That is,
returns imag(get_eigenvalues())
- the returned vector has dim() rows
- the eigenvalues are not sorted in any particular way
!*/
const complex_matrix_type get_v (
) const;
/*!
ensures
- returns the eigenvector matrix V that is
dim() rows by dim() columns
- Each column in V is one of the eigenvectors of the input
matrix
!*/
const complex_matrix_type get_d (
) const;
/*!
ensures
- returns a matrix D such that:
- D.nr() == dim()
- D.nc() == dim()
- diag(D) == get_eigenvalues()
(i.e. the diagonal of D contains all the eigenvalues in the input matrix)
- all off diagonal elements of D are set to 0
!*/
const matrix_type& get_pseudo_v (
) const;
/*!
ensures
- returns a matrix that is dim() rows by dim() columns
- Let A denote the input matrix given to this object's constructor.
- if (A has any imaginary eigenvalues) then
- returns the pseudo-eigenvector matrix V
- The matrix V returned by this function is structured such that:
- A*V == V*get_pseudo_d()
- else
- returns the eigenvector matrix V with A's eigenvectors as
the columns of V
- A*V == V*diagm(get_real_eigenvalues())
!*/
const matrix_type get_pseudo_d (
) const;
/*!
ensures
- The returned matrix is dim() rows by dim() columns
- Computes and returns the block diagonal eigenvalue matrix.
If the original matrix A is not symmetric, then the eigenvalue
matrix D is block diagonal with the real eigenvalues in 1-by-1
blocks and any complex eigenvalues,
a + i*b, in 2-by-2 blocks, (a, b; -b, a). That is, if the complex
eigenvalues look like
u + iv . . . . .
. u - iv . . . .
. . a + ib . . .
. . . a - ib . .
. . . . x .
. . . . . y
Then D looks like
u v . . . .
-v u . . . .
. . a b . .
. . -b a . .
. . . . x .
. . . . . y
This keeps V (The V you get from get_pseudo_v()) a real matrix in both
symmetric and non-symmetric cases, and A*V = V*D.
- the eigenvalues are not sorted in any particular way
!*/
};
// ----------------------------------------------------------------------------------------
}
#endif // DLIB_MATRIx_LA_FUNCTS_ABSTRACT_