GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-09-16 11:23:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

Document Tridiagonalization class, remove unused types.

Browse Source
This commit is contained in:
Jitse Niesen 2010-05-01 17:52:16 +01:00
parent d9c1224133
commit afed0ef90d
7 changed files with 263 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

250
Eigen/src/Eigenvalues/Tridiagonalization.h
View File

@ -30,14 +30,32 @@
*
* \class Tridiagonalization
*
* \brief Trigiagonal decomposition of a selfadjoint matrix
* \brief Tridiagonal decomposition of a selfadjoint matrix
*
* \param MatrixType the type of the matrix of which we are performing the tridiagonalization
* \tparam _MatrixType the type of the matrix of which we are computing the
* tridiagonal decomposition; this is expected to be an instantiation of the
* Matrix class template.
*
* This class performs a tridiagonal decomposition of a selfadjoint matrix \f$ A \f$ such that:
* \f$ A = Q T Q^* \f$ where \f$ Q \f$ is unitary and \f$ T \f$ a real symmetric tridiagonal matrix.
*
* \sa MatrixBase::tridiagonalize()
* A tridiagonal matrix is a matrix which has nonzero elements only on the
* main diagonal and the first diagonal below and above it. The Hessenberg
* decomposition of a selfadjoint matrix is in fact a tridiagonal
* decomposition. This class is used in SelfAdjointEigenSolver to compute the
* eigenvalues and eigenvectors of a selfadjoint matrix.
*
* Call the function compute() to compute the tridiagonal decomposition of a
* given matrix. Alternatively, you can use the Tridiagonalization(const MatrixType&)
* constructor which computes the tridiagonal Schur decomposition at
* construction time. Once the decomposition is computed, you can use the
* matrixQ() and matrixT() functions to retrieve the matrices Q and T in the
* decomposition.
*
* The documentation of Tridiagonalization(const MatrixType&) contains an
* example of the typical use of this class.
*
* \sa class HessenbergDecomposition, class SelfAdjointEigenSolver
*/
template<typename _MatrixType> class Tridiagonalization
{
@ -46,21 +64,18 @@ template<typename _MatrixType> class Tridiagonalization
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar;
typedef typename ei_packet_traits<Scalar>::type Packet;
enum {
Size = MatrixType::RowsAtCompileTime,
SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
Options = MatrixType::Options,
MaxSize = MatrixType::MaxRowsAtCompileTime,
MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1,
PacketSize = ei_packet_traits<Scalar>::size
MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
};
typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
typedef typename ei_plain_col_type<MatrixType, RealScalar>::type DiagonalType;
typedef Matrix<RealScalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> SubDiagonalType;
typedef typename ei_plain_row_type<MatrixType>::type RowVectorType;
typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
typename Diagonal<MatrixType,0>::RealReturnType,
@ -74,22 +89,53 @@ template<typename _MatrixType> class Tridiagonalization
Block<MatrixType,SizeMinusOne,SizeMinusOne>,0 >
>::ret SubDiagonalReturnType;
/** This constructor initializes a Tridiagonalization object for
* further use with Tridiagonalization::compute()
/** \brief Default constructor.
*
* \param [in] size Positive integer, size of the matrix whose tridiagonal
* decomposition will be computed.
*
* The default constructor is useful in cases in which the user intends to
* perform decompositions via compute(). The \p size parameter is only
* used as a hint. It is not an error to give a wrong \p size, but it may
* impair performance.
*
* \sa compute() for an example.
*/
Tridiagonalization(int size = Size==Dynamic ? 2 : Size)
: m_matrix(size,size), m_hCoeffs(size-1)
{}
/** \brief Constructor; computes tridiagonal decomposition of given matrix.
*
* \param[in] matrix Selfadjoint matrix whose tridiagonal decomposition
* is to be computed.
*
* This constructor calls compute() to compute the tridiagonal decomposition.
*
* Example: \include Tridiagonalization_Tridiagonalization_MatrixType.cpp
* Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out
*/
Tridiagonalization(const MatrixType& matrix)
: m_matrix(matrix), m_hCoeffs(matrix.cols()-1)
{
_compute(m_matrix, m_hCoeffs);
}
/** Computes or re-compute the tridiagonalization for the matrix \a matrix.
/** \brief Computes tridiagonal decomposition of given matrix.
*
* \param[in] matrix Selfadjoint matrix whose tridiagonal decomposition
* is to be computed.
*
* This method allows to re-use the allocated data.
* The tridiagonal decomposition is computed by bringing the columns of
* the matrix successively in the required form using Householder
* reflections. The cost is \f$ 4n^3/3 \f$ flops, where \f$ n \f$ denotes
* the size of the given matrix.
*
* This method reuses of the allocated data in the Tridiagonalization
* object, if the size of the matrix does not change.
*
* Example: \include Tridiagonalization_compute.cpp
* Output: \verbinclude Tridiagonalization_compute.out
*/
void compute(const MatrixType& matrix)
{
@ -98,74 +144,191 @@ template<typename _MatrixType> class Tridiagonalization
_compute(m_matrix, m_hCoeffs);
}
/** \returns the householder coefficients allowing to
* reconstruct the matrix Q from the packed data.
/** \brief Returns the Householder coefficients.
*
* \sa packedMatrix()
* \returns a const reference to the vector of Householder coefficients
*
* \pre Either the constructor Tridiagonalization(const MatrixType&) or
* the member function compute(const MatrixType&) has been called before
* to compute the tridiagonal decomposition of a matrix.
*
* The Householder coefficients allow the reconstruction of the matrix
* \f$ Q \f$ in the tridiagonal decomposition from the packed data.
*
* Example: \include Tridiagonalization_householderCoefficients.cpp
* Output: \verbinclude Tridiagonalization_householderCoefficients.out
*
* \sa packedMatrix(), \ref Householder_Module "Householder module"
*/
inline CoeffVectorType householderCoefficients(void) const { return m_hCoeffs; }
inline CoeffVectorType householderCoefficients() const { return m_hCoeffs; }
/** \returns the internal result of the decomposition.
/** \brief Returns the internal representation of the decomposition
*
* \returns a const reference to a matrix with the internal representation
* of the decomposition.
*
* \pre Either the constructor Tridiagonalization(const MatrixType&) or
* the member function compute(const MatrixType&) has been called before
* to compute the tridiagonal decomposition of a matrix.
*
* The returned matrix contains the following information:
* - the strict upper part is equal to the input matrix A
* - the diagonal and lower sub-diagonal represent the tridiagonal symmetric matrix (real).
* - the rest of the lower part contains the Householder vectors that, combined with
* Householder coefficients returned by householderCoefficients(),
* allows to reconstruct the matrix Q as follow:
* Q = H_{N-1} ... H_1 H_0
* where the matrices H are the Householder transformations:
* H_i = (I - h_i * v_i * v_i')
* where h_i == householderCoefficients()[i] and v_i is a Householder vector:
* v_i = [ 0, ..., 0, 1, M(i+2,i), ..., M(N-1,i) ]
* - the strict upper triangular part is equal to the input matrix A.
* - the diagonal and lower sub-diagonal represent the real tridiagonal
* symmetric matrix T.
* - the rest of the lower part contains the Householder vectors that,
* combined with Householder coefficients returned by
* householderCoefficients(), allows to reconstruct the matrix Q as
* \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
* Here, the matrices \f$ H_i \f$ are the Householder transformations
* \f$ H_i = (I - h_i v_i v_i^T) \f$
* where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and
* \f$ v_i \f$ is the Householder vector defined by
* \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$
* with M the matrix returned by this function.
*
* See LAPACK for further details on this packed storage.
*
* Example: \include Tridiagonalization_packedMatrix.cpp
* Output: \verbinclude Tridiagonalization_packedMatrix.out
*
* \sa householderCoefficients()
*/
inline const MatrixType& packedMatrix(void) const { return m_matrix; }
inline const MatrixType& packedMatrix() const { return m_matrix; }
/** \brief Reconstructs the unitary matrix Q in the decomposition
*
* \returns the matrix Q
*
* \pre Either the constructor Tridiagonalization(const MatrixType&) or
* the member function compute(const MatrixType&) has been called before
* to compute the tridiagonal decomposition of a matrix.
*
* This function reconstructs the matrix Q from the Householder
* coefficients and the packed matrix stored internally. This
* reconstruction requires \f$ 4n^3 / 3 \f$ flops.
*
* \sa Tridiagonalization(const MatrixType&) for an example,
* matrixT(), matrixQInPlace()
*/
MatrixType matrixQ() const;
/** \brief Reconstructs the unitary matrix Q in the decomposition
*
* This is an in-place variant of matrixQ() which avoids the copy.
* This function will probably be deleted soon.
*/
template<typename QDerived> void matrixQInPlace(MatrixBase<QDerived>* q) const;
/** \brief Constructs the tridiagonal matrix T in the decomposition
*
* \returns the matrix T
*
* \pre Either the constructor Tridiagonalization(const MatrixType&) or
* the member function compute(const MatrixType&) has been called before
* to compute the tridiagonal decomposition of a matrix.
*
* This function copies the matrix T from internal data. The diagonal and
* subdiagonal of the packed matrix as returned by packedMatrix()
* represents the matrix T. It may sometimes be sufficient to directly use
* the packed matrix or the vector expressions returned by diagonal()
* and subDiagonal() instead of creating a new matrix with this function.
*
* \sa Tridiagonalization(const MatrixType&) for an example,
* matrixQ(), packedMatrix(), diagonal(), subDiagonal()
*/
MatrixType matrixT() const;
const DiagonalReturnType diagonal(void) const;
const SubDiagonalReturnType subDiagonal(void) const;
/** \brief Returns the diagonal of the tridiagonal matrix T in the decomposition.
*
* \returns expression representing the diagonal of T
*
* \pre Either the constructor Tridiagonalization(const MatrixType&) or
* the member function compute(const MatrixType&) has been called before
* to compute the tridiagonal decomposition of a matrix.
*
* Example: \include Tridiagonalization_diagonal.cpp
* Output: \verbinclude Tridiagonalization_diagonal.out
*
* \sa matrixT(), subDiagonal()
*/
const DiagonalReturnType diagonal() const;
/** \brief Returns the subdiagonal of the tridiagonal matrix T in the decomposition.
*
* \returns expression representing the subdiagonal of T
*
* \pre Either the constructor Tridiagonalization(const MatrixType&) or
* the member function compute(const MatrixType&) has been called before
* to compute the tridiagonal decomposition of a matrix.
*
* \sa diagonal() for an example, matrixT()
*/
const SubDiagonalReturnType subDiagonal() const;
/** \brief Performs a full decomposition in place
*
* \param[in,out] mat On input, the selfadjoint matrix whose tridiagonal
* decomposition is to be computed. On output, the orthogonal matrix Q
* in the decomposition if \p extractQ is true.
* \param[out] diag The diagonal of the tridiagonal matrix T in the
* decomposition.
* \param[out] subdiag The subdiagonal of the tridiagonal matrix T in
* the decomposition.
* \param[in] extractQ If true, the orthogonal matrix Q in the
* decomposition is computed and stored in \p mat.
*
* Compute the tridiagonal matrix of \p mat in place. The tridiagonal
* matrix T is passed to the output parameters \p diag and \p subdiag. If
* \p extractQ is true, then the orthogonal matrix Q is passed to \p mat.
*
* The vectors \p diag and \p subdiag are not resized. The function
* assumes that they are already of the correct size. The length of the
* vector \p diag should equal the number of rows in \p mat, and the
* length of the vector \p subdiag should be one left.
*
* This implementation contains an optimized path for real 3-by-3 matrices
* which is especially useful for plane fitting.
*
* \note Notwithstanding the name, the current implementation copies
* \p mat to a temporary matrix and uses that matrix to compute the
* decomposition.
*
* Example (this uses the same matrix as the example in
* Tridiagonalization(const MatrixType&)):
* \include Tridiagonalization_decomposeInPlace.cpp
* Output: \verbinclude Tridiagonalization_decomposeInPlace.out
*
* \sa Tridiagonalization(const MatrixType&), compute()
*/
static void decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ = true);
static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs);
protected:
static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs);
static void _decomposeInPlace3x3(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ = true);
MatrixType m_matrix;
CoeffVectorType m_hCoeffs;
};
/** \returns an expression of the diagonal vector */
template<typename MatrixType>
const typename Tridiagonalization<MatrixType>::DiagonalReturnType
Tridiagonalization<MatrixType>::diagonal(void) const
Tridiagonalization<MatrixType>::diagonal() const
{
return m_matrix.diagonal();
}
/** \returns an expression of the sub-diagonal vector */
template<typename MatrixType>
const typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
Tridiagonalization<MatrixType>::subDiagonal(void) const
Tridiagonalization<MatrixType>::subDiagonal() const
{
int n = m_matrix.rows();
return Block<MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal();
}
/** constructs and returns the tridiagonal matrix T.
* Note that the matrix T is equivalent to the diagonal and sub-diagonal of the packed matrix.
* Therefore, it might be often sufficient to directly use the packed matrix, or the vector
* expressions returned by diagonal() and subDiagonal() instead of creating a new matrix.
*/
template<typename MatrixType>
typename Tridiagonalization<MatrixType>::MatrixType
Tridiagonalization<MatrixType>::matrixT(void) const
Tridiagonalization<MatrixType>::matrixT() const
{
// FIXME should this function (and other similar ones) rather take a matrix as argument
// and fill it ? (to avoid temporaries)
@ -223,10 +386,9 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
}
}
/** reconstructs and returns the matrix Q */
template<typename MatrixType>
typename Tridiagonalization<MatrixType>::MatrixType
Tridiagonalization<MatrixType>::matrixQ(void) const
Tridiagonalization<MatrixType>::matrixQ() const
{
MatrixType matQ;
matrixQInPlace(&matQ);
@ -240,6 +402,7 @@ void Tridiagonalization<MatrixType>::matrixQInPlace(MatrixBase<QDerived>* q) con
QDerived& matQ = q->derived();
int n = m_matrix.rows();
matQ = MatrixType::Identity(n,n);
typedef typename ei_plain_row_type<MatrixType>::type RowVectorType;
RowVectorType aux(n);
for (int i = n-2; i>=0; i--)
{
@ -248,7 +411,6 @@ void Tridiagonalization<MatrixType>::matrixQInPlace(MatrixBase<QDerived>* q) con
}
}
/** Performs a full decomposition in place */
template<typename MatrixType>
void Tridiagonalization<MatrixType>::decomposeInPlace(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ)
{

11
doc/snippets/Tridiagonalization_Tridiagonalization_MatrixType.cpp Normal file
View File

@ -0,0 +1,11 @@
MatrixXd X = MatrixXd::Random(5,5);
MatrixXd A = X + X.transpose();
cout << "Here is a random symmetric 5x5 matrix:" << endl << A << endl << endl;
Tridiagonalization<MatrixXd> triOfA(A);
cout << "The orthogonal matrix Q is:" << endl << triOfA.matrixQ() << endl;
cout << "The tridiagonal matrix T is:" << endl << triOfA.matrixT() << endl << endl;
MatrixXd Q = triOfA.matrixQ();
MatrixXd T = triOfA.matrixT();
cout << "Q * T * Q^T = " << endl << Q * T * Q.transpose() << endl;

9
doc/snippets/Tridiagonalization_compute.cpp Normal file
View File

@ -0,0 +1,9 @@
Tridiagonalization<MatrixXf> tri;
MatrixXf X = MatrixXf::Random(4,4);
MatrixXf A = X + X.transpose();
tri.compute(A);
cout << "The matrix T in the tridiagonal decomposition of A is: " << endl;
cout << tri.matrixT() << endl;
tri.compute(2*A); // re-use tri to compute eigenvalues of 2A
cout << "The matrix T in the tridiagonal decomposition of 2A is: " << endl;
cout << tri.matrixT() << endl;

10
doc/snippets/Tridiagonalization_decomposeInPlace.cpp Normal file
View File

@ -0,0 +1,10 @@
MatrixXd X = MatrixXd::Random(5,5);
MatrixXd A = X + X.transpose();
cout << "Here is a random symmetric 5x5 matrix:" << endl << A << endl << endl;
VectorXd diag(5);
VectorXd subdiag(4);
Tridiagonalization<MatrixXd>::decomposeInPlace(A, diag, subdiag);
cout << "The orthogonal matrix Q is:" << endl << A << endl;
cout << "The diagonal of the tridiagonal matrix T is:" << endl << diag << endl;
cout << "The subdiagonal of the tridiagonal matrix T is:" << endl << subdiag << endl;

13
doc/snippets/Tridiagonalization_diagonal.cpp Normal file
View File

@ -0,0 +1,13 @@
MatrixXcd X = MatrixXcd::Random(4,4);
MatrixXcd A = X + X.adjoint();
cout << "Here is a random self-adjoint 4x4 matrix:" << endl << A << endl << endl;
Tridiagonalization<MatrixXcd> triOfA(A);
MatrixXcd T = triOfA.matrixT();
cout << "The tridiagonal matrix T is:" << endl << triOfA.matrixT() << endl << endl;
cout << "We can also extract the diagonals of T directly ..." << endl;
VectorXd diag = triOfA.diagonal();
cout << "The diagonal is:" << endl << diag << endl;
VectorXd subdiag = triOfA.subDiagonal();
cout << "The subdiagonal is:" << endl << subdiag << endl;

6
doc/snippets/Tridiagonalization_householderCoefficients.cpp Normal file
View File

@ -0,0 +1,6 @@
Matrix4d X = Matrix4d::Random(4,4);
Matrix4d A = X + X.transpose();
cout << "Here is a random symmetric 4x4 matrix:" << endl << A << endl;
Tridiagonalization<Matrix4d> triOfA(A);
Vector3d hc = triOfA.householderCoefficients();
cout << "The vector of Householder coefficients is:" << endl << hc << endl;

8
doc/snippets/Tridiagonalization_packedMatrix.cpp Normal file
View File

@ -0,0 +1,8 @@
Matrix4d X = Matrix4d::Random(4,4);
Matrix4d A = X + X.transpose();
cout << "Here is a random symmetric 4x4 matrix:" << endl << A << endl;
Tridiagonalization<Matrix4d> triOfA(A);
Matrix4d pm = triOfA.packedMatrix();
cout << "The packed matrix M is:" << endl << pm << endl;
cout << "The diagonal and subdiagonal corresponds to the matrix T, which is:"
<< endl << triOfA.matrixT() << endl;