coco
1006 lines
38 KiB
C++
1006 lines
38 KiB
C++
// 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 sparse_vector_type,
|
|
typename T
|
|
>
|
|
void svd_fast (
|
|
const std::vector<sparse_vector_type>& A,
|
|
matrix<T>& w,
|
|
matrix<T>& v,
|
|
unsigned long l,
|
|
unsigned long q = 1
|
|
);
|
|
/*!
|
|
This function is identical to the above svd_fast() except it doesn't compute u.
|
|
!*/
|
|
|
|
// ----------------------------------------------------------------------------------------
|
|
|
|
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_
|
|
|