diff options
Diffstat (limited to 'ml/dlib/dlib/matrix/matrix_la_abstract.h')
| -rw-r--r-- | ml/dlib/dlib/matrix/matrix_la_abstract.h | 1005 |
1 files changed, 1005 insertions, 0 deletions
diff --git a/ml/dlib/dlib/matrix/matrix_la_abstract.h b/ml/dlib/dlib/matrix/matrix_la_abstract.h new file mode 100644 index 000000000..df6a5fd33 --- /dev/null +++ b/ml/dlib/dlib/matrix/matrix_la_abstract.h @@ -0,0 +1,1005 @@ +// 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_ + |