Interfaces to Lapack and BLAS for Eigen-like arrays.

Functions
template<typename eigenT1, typename eigenT2>
void	mx::math::eigenSYRK (eigenT1 &cv, const eigenT2 &ims)
	Calculates the lower triangular part of the covariance matrix of ims.
template<typename arrT>
MXLAPACK_INT	mx::math::eigenSYEVR (arrT &eigvec, arrT &eigval, arrT &X, int ev0=0, int ev1=-1, char UPLO='L', syevrMem< typename arrT::Scalar > *mem=0)
	Calculate select eigenvalues and eigenvectors of an Eigen Array.
template<typename _evCalcT = double, typename eigenT>
MXLAPACK_INT	mx::math::calcEigenVecs (eigenT &evecs, eigenT &evals, eigenT &cv, int nVecs=0, bool normalize=false, bool check=false, syevrMem< _evCalcT > mem=0, double t_eigenv=nullptr)
	Calculate the eigenvectors and eigenvalues given a triangular matrix.
template<typename _evCalcT = double, typename eigenT, typename eigenT1>
MXLAPACK_INT	mx::math::calcKLModes (eigenT &klModes, eigenT &cv, const eigenT1 &Rims, int n_modes=0, syevrMem< _evCalcT > mem=0, double t_eigenv=nullptr, double *t_klim=nullptr)
	Calculate the K-L modes, or principle components, given a covariance matrix.
template<typename dataT>
MXLAPACK_INT	mx::math::eigenGESDD (Eigen::Array< dataT, -1, -1 > &U, Eigen::Array< dataT, -1, -1 > &S, Eigen::Array< dataT, -1, -1 > &VT, Eigen::Array< dataT, -1, -1 > &A)
	Compute the SVD of an Eigen::Array using LAPACK's xgesdd.
template<typename dataT>
int	mx::math::eigenPseudoInverse (Eigen::Array< dataT, -1, -1 > &PInv, dataT &condition, int &nRejected, Eigen::Array< dataT, -1, -1 > &U, Eigen::Array< dataT, -1, -1 > &S, Eigen::Array< dataT, -1, -1 > &VT, int minMN, dataT &maxCondition, dataT alpha=0, int interact=MX_PINV_NO_INTERACT)
	Calculate the pseudo-inverse of a patrix given its SVD.
template<typename dataT>
int	mx::math::eigenPseudoInverse (Eigen::Array< dataT, -1, -1 > &PInv, dataT &condition, int &nRejected, Eigen::Array< dataT, -1, -1 > &U, Eigen::Array< dataT, -1, -1 > &S, Eigen::Array< dataT, -1, -1 > &VT, Eigen::Array< dataT, -1, -1 > &A, dataT &maxCondition, dataT alpha=0, int interact=MX_PINV_NO_INTERACT)
	Calculate the pseudo-inverse of a patrix using the SVD.
template<typename dataT>
int	mx::math::eigenPseudoInverse (Eigen::Array< dataT, -1, -1 > &PInv, dataT &condition, int &nRejected, Eigen::Array< dataT, -1, -1 > &A, dataT &maxCondition, dataT alpha=0, int interact=MX_PINV_NO_INTERACT)
	Calculate the pseudo-inverse of a matrix using the SVD.

Function Documentation

◆ calcEigenVecs()

template<typename _evCalcT = double, typename eigenT>

MXLAPACK_INT mx::math::calcEigenVecs	(	eigenT &	evecs,
		eigenT &	evals,
		eigenT &	cv,
		int	nVecs = 0,
		bool	normalize = false,
		bool	check = false,
		syevrMem< _evCalcT > *	mem = 0,
		double *	t_eigenv = nullptr )

Calculate the eigenvectors and eigenvalues given a triangular matrix.

Eigen-decomposition of the matrix is performed using eigenSYEVR().

Template Parameters

evCalcT	is the type in which to perform eigen-decomposition, which may be different from the input array and output arrays (which must be the same).
eigenT	is a 2D Eigen-like type

Parameters

[out]	evecs	on exit contains the eigen vectors
[out]	evals	on exit contains the eigen vectors
[in]	cv	a lower-triangle (in the Lapack sense) square covariance matrix.
[in]	nVecs	[opt] The maximum number of modes to solve for. If 0 all modes are solved for.
[in]	normalize	[opt] flag specifying whether or not to normalize the eigenvectors.
[in]	check	[opt] flag specifying whether or not to check the eigenvalues/vectors for validity. Requires normalize=true.
[in]	mem	[opt] A memory structure which can be re-used by SYEVR for efficiency.
[out]	t_eigenv	[opt] if not null, will be filled in with the time taken to calculate eigenvalues.

Definition at line 537 of file eigenLapack.hpp.

References eigenSYEVR(), mx::sys::get_curr_time(), and isFinite().

Referenced by calcKLModes(), unitTest::math_eigenLapack_test::TEST_CASE(), unitTest::math_eigenLapack_test::TEST_CASE(), and unitTest::math_eigenLapack_test::TEST_CASE().

◆ calcKLModes()

template<typename _evCalcT = double, typename eigenT, typename eigenT1>

MXLAPACK_INT mx::math::calcKLModes	(	eigenT &	klModes,
		eigenT &	cv,
		const eigenT1 &	Rims,
		int	n_modes = 0,
		syevrMem< _evCalcT > *	mem = 0,
		double *	t_eigenv = nullptr,
		double *	t_klim = nullptr )

Calculate the K-L modes, or principle components, given a covariance matrix.

Eigen-decomposition of the covariance matrix is performed using eigenSYEVR().

Template Parameters

evCalcT	is the type in which to perform eigen-decomposition.
eigenT	is a 2D Eigen-like type
eigenT1	is a 2D Eigen-like type.

Parameters

[out]	klModes	on exit contains the K-L modes (or P.C.s)
[in]	cv	a lower-triangle (in the Lapack sense) square covariance matrix.
[in]	Rims	The reference data. cv.rows() == Rims.cols().
[in]	n_modes	[opt] Tbe maximum number of modes to solve for. If 0 all modes are solved for.
[in]	mem	[opt] A memory structure which can be re-used by SYEVR for efficiency.
[out]	t_eigenv	[opt] if not null, will be filled in with the time taken to calculate eigenvalues.
[out]	t_klim	[opt] if not null, will be filled in with the time taken to calculate the KL modes.

Definition at line 679 of file eigenLapack.hpp.

References calcEigenVecs(), gemm(), and mx::sys::get_curr_time().

Referenced by unitTest::math_eigenLapack_test::TEST_CASE(), and unitTest::math_eigenLapack_test::TEST_CASE().

◆ eigenGESDD()

template<typename dataT>

MXLAPACK_INT mx::math::eigenGESDD	(	Eigen::Array< dataT, -1, -1 > &	U,
		Eigen::Array< dataT, -1, -1 > &	S,
		Eigen::Array< dataT, -1, -1 > &	VT,
		Eigen::Array< dataT, -1, -1 > &	A )

Compute the SVD of an Eigen::Array using LAPACK's xgesdd.

Computes the SVD of A, \( A = U S V^T \).

Returns: 0 on success; -i on error in ith parameter (from LAPACK xgesdd); >0 did not converge (from LAPACK xgesdd)

Template Parameters

dataT is either float or double.

Parameters

[out]	U	the A.rows() x A.rows() left matrix
[out]	S	the A.cols() x 1 matrix of singular values
[out]	VT	the A.cols() x A.cols() right matrix, note this is the transpose.
[in]	A	the input matrix to be decomposed

Definition at line 794 of file eigenLapack.hpp.

References gesdd().

Referenced by eigenPseudoInverse().

◆ eigenPseudoInverse() [1/3]

template<typename dataT>

int mx::math::eigenPseudoInverse	(	Eigen::Array< dataT, -1, -1 > &	PInv,
		dataT &	condition,
		int &	nRejected,
		Eigen::Array< dataT, -1, -1 > &	A,
		dataT &	maxCondition,
		dataT	alpha = 0,
		int	interact = MX_PINV_NO_INTERACT )

Calculate the pseudo-inverse of a matrix using the SVD.

First computes the SVD of A, \( A = U S V^T \), using eigenGESDD. Then the psuedo-inverse is calculated as \( A^+ = V S^+ U^T\). This interface does not provide access to U, S and VT.

The parameter interact is intepreted as a bitmask. The values can be

MX_PINV_PLOT, which will cause a plot to be displayed of the singular values
MX_PINV_ASK, which will ask the user for a max. condition number using stdin
MX_PINV_ASK_NMODES, which will ask the user for a max number of modes to include using stdin. Overrides MX_PINV_ASK. If interact is 0 then no interaction is used and maxCondition controls the inversion.
- Template Parameters
  
  dataT is either float or double.
  
  This is an overloaded member function, provided for convenience. It differs from the above function only in what argument(s) it accepts.

Parameters

[out]	PInv	The pseudo-inverse of A
[out]	condition	The final condition number.
[out]	nRejected	The number of eigenvectors rejected
[in]	A	The matrix to invert, will be altered!
[in]	maxCondition	If > 0, the maximum condition number desired. If <0 the number of modes to keep. Used to threshold the singular values. Set to 0 to include all eigenvalues/vectors. Ignored if interactive.
[in]	alpha	[opt] the Tikhonov regularization value, as a fraction of the highest singular value. If alpha < 0, then it is treated as a (positive) floor (as a fraction of highest singular value) for the singular values, which is not the same as Tikhonov (alpha > 0).
[in]	interact	[opt] a bitmask controlling interaction. See above.

Definition at line 1102 of file eigenLapack.hpp.

References eigenPseudoInverse().

◆ eigenPseudoInverse() [2/3]

template<typename dataT>

int mx::math::eigenPseudoInverse	(	Eigen::Array< dataT, -1, -1 > &	PInv,
		dataT &	condition,
		int &	nRejected,
		Eigen::Array< dataT, -1, -1 > &	U,
		Eigen::Array< dataT, -1, -1 > &	S,
		Eigen::Array< dataT, -1, -1 > &	VT,
		Eigen::Array< dataT, -1, -1 > &	A,
		dataT &	maxCondition,
		dataT	alpha = 0,
		int	interact = MX_PINV_NO_INTERACT )

Calculate the pseudo-inverse of a patrix using the SVD.

First computes the SVD of A, \( A = U S V^T \), using eigenGESDD. Then the psuedo-inverse is calculated as \( A^+ = V S^+ U^T\).

The parameter interact is intepreted as a bitmask. The values can be

MX_PINV_PLOT, which will cause a plot to be displayed of the singular values
MX_PINV_ASK, which will ask the user for a max. condition number using stdin
MX_PINV_ASK_NMODES, which will ask the user for a max number of modes to include using stdin. Overrides MX_PINV_ASK. If interact is 0 then no interaction is used and maxCondition controls the inversion.

Template Parameters

dataT is either float or double.

Parameters

[out]	PInv	The pseudo-inverse of A
[out]	condition	The final condition number.
[out]	nRejected	The number of eigenvectors rejected
[out]	U	the A.rows() x A.rows() left matrix
[out]	S	the A.cols() x 1 matrix of singular values
[out]	VT	the A.cols() x A.cols() right matrix, note this is the transpose.
[in]	A	The matrix to invert. This will be modified!
[in]	maxCondition	If > 0, the maximum condition number desired. If <0 the number of modes to keep. Used to threshold the singular values. Set to 0 to include all eigenvalues/vectors. Ignored if interactive.
[in]	alpha	[opt] the Tikhonov regularization value, as a fraction of the highest singular value. If alpha < 0, then it is treated as a (positive) floor (as a fraction of highest singular value) for the singular values, which is not the same as Tikhonov (alpha > 0).
[in]	interact	[opt] a bitmask controlling interaction. See above.

Definition at line 1046 of file eigenLapack.hpp.

References eigenGESDD(), and eigenPseudoInverse().

◆ eigenPseudoInverse() [3/3]

template<typename dataT>

int mx::math::eigenPseudoInverse	(	Eigen::Array< dataT, -1, -1 > &	PInv,
		dataT &	condition,
		int &	nRejected,
		Eigen::Array< dataT, -1, -1 > &	U,
		Eigen::Array< dataT, -1, -1 > &	S,
		Eigen::Array< dataT, -1, -1 > &	VT,
		int	minMN,
		dataT &	maxCondition,
		dataT	alpha = 0,
		int	interact = MX_PINV_NO_INTERACT )

Calculate the pseudo-inverse of a patrix given its SVD.

Given the SVD of A, \( A = U S V^T \), as calculated by eigenGESDD the psuedo-inverse is calculated as \( A^+ = V S^+ U^T\).

The parameter interact is intepreted as a bitmask. The values can be

MX_PINV_PLOT, which will cause a plot to be displayed of the singular values
MX_PINV_ASK, which will ask the user for a max. condition number using stdin
MX_PINV_ASK_NMODES, which will ask the user for a max number of modes to include using stdin. Overrides MX_PINV_ASK. If interact is 0 then no interaction is used and maxCondition controls the inversion.

Template Parameters

dataT is either float or double.

Parameters

[out]	PInv	The pseudo-inverse of A
[out]	condition	The final condition number.
[out]	nRejected	The number of eigenvectors rejected
[in]	U	the A.rows() x A.rows() left matrix
[in]	S	the A.cols() x 1 matrix of singular values
[in]	VT	the A.cols() x A.cols() right matrix, note this is the transpose.
[in]	minMN	The minimum size of the matrix to invert.
[in]	maxCondition	If > 0, the maximum condition number desired. If <0 the number of modes to keep. Used to threshold the singular values. Set to 0 to include all eigenvalues/vectors. Ignored if interactive.
[in]	alpha	[opt] the Tikhonov regularization value, as a fraction of the highest singular value. If alpha < 0, then it is treated as a (positive) floor (as a fraction of highest singular value) for the singular values, which is not the same as Tikhonov (alpha > 0).
[in]	interact	[opt] a bitmask controlling interaction. See above.

Definition at line 854 of file eigenLapack.hpp.

References mx::math::gnuPlot::command(), mx::math::gnuPlot::logy(), and mx::math::gnuPlot::plot().

Referenced by eigenPseudoInverse(), eigenPseudoInverse(), and mx::AO::ifPInv().

◆ eigenSYEVR()

template<typename arrT>

MXLAPACK_INT mx::math::eigenSYEVR	(	arrT &	eigvec,
		arrT &	eigval,
		arrT &	X,
		int	ev0 = 0,
		int	ev1 = -1,
		char	UPLO = 'L',
		syevrMem< typename arrT::Scalar > *	mem = 0 )

Calculate select eigenvalues and eigenvectors of an Eigen Array.

Uses the templateLapack wrapper for syevr.

Template Parameters

arrT	is the eigen-like type containing the data

Returns: -1 for invalid matrix/range geometry or invalid workspace-query results.; -1000 on an malloc allocation error.; the return code from syevr (info) otherwise.

Parameters

[out]	eigvec	will contain the eigenvectors as columns
[out]	eigval	will contain the eigenvalues
[in]	X	is a square matrix which is either upper or lower (default) triangular
[in]	ev0	[opt] is the first desired eigenvalue (default = 0)
[in]	ev1	[opt] exclusive upper bound of the desired eigenvalue indices. If -1 all eigenvalues are returned.
[in]	UPLO	[opt] specifies whether X is upper ('U') or lower ('L') triangular. Default is ('L').
[in]	mem	[opt] holds the working memory arrays, can be re-passed to avoid unnecessary re-allocations

Definition at line 360 of file eigenLapack.hpp.

References mx::allocerr, isFinite(), lamch(), mx::lapackerr, and mx::internal::mxlib_error_report().

Referenced by calcEigenVecs(), unitTest::math_eigenLapack_test::TEST_CASE(), unitTest::math_eigenLapack_test::TEST_CASE(), unitTest::math_eigenLapack_test::TEST_CASE(), and unitTest::math_eigenLapack_test::TEST_CASE().

◆ eigenSYRK()

template<typename eigenT1, typename eigenT2>

void mx::math::eigenSYRK	(	eigenT1 &	cv,
		const eigenT2 &	ims )

Calculates the lower triangular part of the covariance matrix of ims.

Uses cblas_ssyrk. cv is resized to ims.cols() X ims.cols(). Calculates \( cv = A^T*A \).

Template Parameters

eigenT1	is the eigen matrix/array type of cv.
eigenT2	is the eigen matrix/array type of ims

Parameters

[out]	cv	is the eigen matrix/array where to store the result
[in]	ims	is the eigen matrix/array (images as columns) to calculate the covariance of

Definition at line 230 of file eigenLapack.hpp.

References syrk().

Referenced by unitTest::math_eigenLapack_test::TEST_CASE().

Functions

Function Documentation

◆ calcEigenVecs()

◆ calcKLModes()

◆ eigenGESDD()

◆ eigenPseudoInverse() [1/3]

◆ eigenPseudoInverse() [2/3]

◆ eigenPseudoInverse() [3/3]

◆ eigenSYEVR()

◆ eigenSYRK()