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