28 #include <Base/Debug/Exception.hh>
29 #include <Base/Common/CompareFloatingPoint.hh>
36 #ifdef BIAS_HAVE_OPENCV
44 bool AbsoluteZeroThreshold) {
46 AbsoluteZeroThreshold_ = AbsoluteZeroThreshold;
47 Compute( M, ZeroThreshold);
50 #ifdef BIAS_HAVE_OPENCV
52 ZeroThreshold_ = ZeroThreshold;
56 CvMat* mat = cvCreateMat(rows,cols,CV_32FC1);
59 for(
unsigned i=0;i<rows;i++){
60 for(
unsigned j=0;j<cols;j++){
61 CV_MAT_ELEM( *mat,
float, i, j ) = (float)M[i][j];
65 CvMat* U = cvCreateMat(rows,rows,CV_32FC1); cvZero(U);
66 CvMat* D = cvCreateMat(rows,cols,CV_32FC1); cvZero(D);
67 CvMat* V = cvCreateMat(cols,cols,CV_32FC1); cvZero(V);
68 cvSVD(mat, D, U, V, CV_SVD_V_T);
70 U_.newsize(rows,rows);
71 VT_.newsize(cols,cols);
73 for(
unsigned i=0;i<rows;i++){
74 for(
unsigned j=0;j<rows;j++){
75 U_[i][j] = (double)CV_MAT_ELEM(*U,
float,i,j);
77 if(i< rows && i < cols){
79 S_[i] = (double)CV_MAT_ELEM(*D,
float,i,i);
82 for(
unsigned i=0;i<cols;i++){
83 for(
unsigned j=0;j<cols;j++){
84 VT_[i][j] = (double)CV_MAT_ELEM(*V,
float,i,j);
94 if (S_[0] < numeric_limits<double>::epsilon())
105 ZeroThreshold_ = ZeroThreshold;
106 return General_Eigenproblem_GeneralMatrix_Lapack( M );
121 if ( S_.size() != 0 && U_.num_rows() != 0 && U_.num_cols() != 0 &&
122 VT_.num_rows() != 0 && VT_.num_cols() != 0 && res == 0) {
124 if (S_[0] < numeric_limits<double>::epsilon()) {
126 }
else Decomposed_=
true;
138 BEXCEPTION(
"Error, there is a bug in this function, use Solve(Matrix, "
139 "vector, Vector instead");
146 if (y.
size() <= U_.num_rows()){
149 if (U_.num_rows() < U_.num_cols()) {
151 for (
int k=0; k< y.
size();k++)
160 for (; i < UTy.
size(); i++)
161 x[i]=(fabs(S_[i]) >= ZeroThreshold_) ? UTy[i]/S_[i] : 0;
163 for (;i<x.size();i++) x[i]=0;
165 return VT_.Transpose() * x;
168 BIASERR(
"size of right hand side of Ax = y is too big. Length(y) = "
174 BIASERR(
"SVD holds no decomposed matrix.");
180 #ifdef BIAS_HAVE_OPENCV
187 CvMat* a = cvCreateMat(rows,cols,CV_32FC1);
190 for(
unsigned i=0;i<rows;i++){
191 for(
unsigned j=0;j<cols;j++){
192 CV_MAT_ELEM( *a,
float, i, j ) = (float)A[i][j];
195 CvMat* b = cvCreateMat(rows,cols,CV_32FC1);
199 double ret = cvInvert(a,b,CV_SVD);
200 if(ret==0.0) BIASWARN(
"Could not invert Matrix");
202 for(
unsigned i=0;i<rows;i++){
203 for(
unsigned j=0;j<cols;j++){
204 B[i][j] = (double)CV_MAT_ELEM( *b,
float, i, j );
226 if (AbsoluteZeroThreshold_) {
227 for (
register int i = 0; i < S_.size(); i++){
228 SdT[i][i] = (S_[i] < ZeroThreshold_ && S_[i] > -ZeroThreshold_)
229 ? (0.0) : (1.0/S_[i]);
232 for (
register int i = 0; i < S_.size(); i++){
233 SdT[i][i]=(S_[i]/S_[0] < ZeroThreshold_ && S_[i]/S_[0] > -ZeroThreshold_)
234 ? (0.0) : (1.0/S_[i]);
255 for (
register int i = 0; i < rank; i++){
256 SdT[i][i] = (1.0/S_[i]);
258 for (
register int i = rank; i < S_.size(); i++){
281 bool valid = (A_.GetCols() == A_.GetRows());
283 const int num = A_.GetCols();
284 for (
int r=0; r<num; r++){
285 for (
int c=r; c<num; c++){
287 valid = valid &&
Equal(A_[r][c], A_[c][r], 1e-10);
289 cout << setprecision(20) << A_[r][c] <<
"\t"<< setprecision(20) << A_[c][r]<<endl;
298 BEXCEPTION(
"matrix is not symmetric ");
302 for (
register int i=0; i<eig.
size(); i++){
305 cout <<
"eigenvalue "<<setprecision(20)<<eig[i]<<endl;
311 BEXCEPTION(
"matrix is not positive definit ");
316 if ((res = Compute(A_))!=0) {
317 BIASERR(
"SVD failed with result "<<res);
324 if (AbsoluteZeroThreshold_) {
325 for (
register int i=0; i<S_.size(); i++){
326 SqrtS[i][i] = (fabs(S_[i])<ZeroThreshold_)?(0.0):(sqrt(S_[i]));
329 for (
register int i=0; i<S_.size(); i++){
330 SqrtS[i][i] = (fabs(S_[i]/S_[0])<ZeroThreshold_)?(0.0):(sqrt(S_[i]));
334 return U_ * SqrtS * VT_;
349 bool valid = (A_.GetCols() == A_.GetRows());
351 const int num = A_.GetCols();
352 for (
int r=0; r<num; r++){
353 for (
int c=r; c<num; c++){
354 valid = valid &&
Equal(A_[r][c], A_[c][r], 1e-12);
356 cout<<
" Invalid matrix ("<<num<<
"x"<<num<<
") element at "
358 cout << setprecision(20) << A_[r][c] <<
"\t"<< setprecision(20)
363 if (!valid) {
break; }
367 BEXCEPTION(
"matrix is not symmetric "<<A_<<endl
368 <<
"Fix this by calling A_.MakeSymmetric() in CALLING class");
373 for (
register int i=0; i<eig.
size(); i++){
387 if (Compute(A_)!=0) {
394 if (AbsoluteZeroThreshold_) {
395 for (
register int i=0; i<S_.size(); i++){
396 SqrtS[i][i] = (fabs(S_[i])<ZeroThreshold_)?(0.0):(sqrt(S_[i]));
399 for (
register int i=0; i<S_.size(); i++){
400 SqrtS[i][i] = (fabs(S_[i]/S_[0])<ZeroThreshold_)?(0.0):(sqrt(S_[i]));
418 int index = S_.size()-1;
421 while (index >= 0 && fabs(S_[index]) < ZeroThreshold_ ) {
427 if (A_.num_cols()>A_.num_rows()){
428 dim+=A_.num_cols()-A_.num_rows();
436 int index = S_.size()-1;
439 const double RefSingValue(S_[0]);
440 if (RefSingValue<=DBL_EPSILON) {
445 while ((S_[index]/RefSingValue)<ZeroThreshold_){
447 if (--index<1)
break;
451 if (A_.num_cols()>A_.num_rows()){
452 dim += A_.num_cols()-A_.num_rows();
460 return AbsNullspaceDim();
464 return RelNullspaceDim();
470 int index = S_.size()-1;
473 while (index >= 0 && fabs(S_[index]) < ZeroThreshold_ ) {
479 if (A_.num_rows()>A_.num_cols()){
480 dim += A_.num_rows()-A_.num_cols();
488 int index = S_.size()-1;
490 if (GetSingularValue(0)==0){
494 while (index>=0 && fabs(S_[index]/maxs)<ZeroThreshold_){
499 if (A_.num_rows()>A_.num_cols()){
500 dim+=A_.num_rows()-A_.num_cols();
508 BIASDOUT(D_SVD_RANK,
"S_.Size(): "<<S_.Size()<<
" NullspaceDim(): "
509 << NullspaceDim()<<
" A_.num_cols(): "<<A_.num_cols()
510 <<
" A_.num_rows(): "<<A_.num_rows()<<
" S_: "<<S_);
511 unsigned int lowdim=A_.num_rows()<A_.num_cols()?A_.num_rows():A_.num_cols();
512 return (lowdim-NullspaceDim());
528 BIASCDOUT(D_SVD_SOLVE,
"S: "<<S<<endl);
540 BIASCDOUT(D_SVD_SOLVE,
"d: "<<d<<endl);
541 if (AbsoluteZeroThreshold_) {
542 for (
int i=0; i<A.
num_cols(); i++) {
543 if (fabs(S[i])>ZeroThreshold_) {
551 for (
int i=0; i<A.
num_cols(); i++) {
552 if (fabs(S[i]/S_[0])>ZeroThreshold_) {
562 BIASDOUT(D_SVD_SOLVE,
"Error in SVD, degenerate case in Solve ("
563 << -res <<
" dim. solution space), S="
575 std::vector<double> eigValues(VT_.GetRows());
579 for (
unsigned int i=0; i<VT_.GetRows(); i++) {
581 eigenvector=VT_.GetRow(i);
582 eigenimage = A_ * eigenvector;
584 eigValues[i] = eigenimage * eigenvector;
586 std::sort(eigValues.begin(),eigValues.end());
int General_Eigenproblem_GeneralMatrix_Lapack(const Matrix< double > &M)
solve the general (non-special) eigenvalue/eigenvector problem of a general (non-symmetric) matrix M ...
const int AbsNullspaceDim() const
returns dim of nullspace using the absolute threshold criterion
const int AbsRightNullspaceDim() const
returns the dim of the matrix's right nullspace, using absolute threshold criterion ...
Subscript num_cols() const
Vector< double > Solve(const Vector< double > &y) const
Matrix< T > Transpose() const
transpose function, storing data destination matrix
unsigned int Rank()
returns the rank of A_
Matrix< double > SqrtT()
returns the square root of a symmetric positive definite matrix M A = sqrt(M) = U * sqrt(S)...
BIAS::Vector< double > GetEigenValues()
Call this after Compute(), returns the eigenvalues of the matrix in ascending order (smallest first)...
int Compute(const Matrix< double > &M, double ZeroThreshold=DEFAULT_DOUBLE_ZERO_THRESHOLD)
set a new matrix and compute its decomposition.
const int RelNullspaceDim() const
compare singular values against greatest, not absolute
void SetZero()
Sets all values to zero.
unsigned int GetRows() const
bool GreaterEqual(const T left, const T right, const T eps=std::numeric_limits< T >::epsilon())
comparison function for floating point values
const int RelLeftNullspaceDim() const
Matrix< double > InvertOpenCV(const Matrix< double > &A)
Calculates new inverse with OpenCV cvInvert.
SVD(const Matrix< double > &M, double ZeroThreshold=DEFAULT_DOUBLE_ZERO_THRESHOLD, bool AbsoluteZeroThreshold=true)
solve the general eigenproblem of the rectangular matrix M or with other words, make the U*S*V^T deco...
int General_singular_value_decomposition(char jobu, char jobvt, const BIAS::Matrix< double > &A, BIAS::Vector< double > &ret_S, BIAS::Matrix< double > &ret_U, BIAS::Matrix< double > &ret_VT)
solve general eigenproblem.
bool Equal(const T left, const T right, const T eps)
comparison function for floating point values See http://www.boost.org/libs/test/doc/components/test_...
const int AbsLeftNullspaceDim() const
returns the dim of the matrix's left nullspace, using absolute threshold criterion ...
unsigned int GetCols() const
int ComputeOpenCV(const Matrix< double > &M, double ZeroThreshold=DEFAULT_DOUBLE_ZERO_THRESHOLD)
Use OpenCV for decomposition as Lapack frequently leads to crashes under windows. ...
Matrix< double > Invert()
returns pseudoinverse of A = U * S * V^T A^+ = V * S^+ * U^T
Matrix< double > Sqrt()
returns the square root of a symmetric positive definite matrix M A = sqrt(M) = U * sqrt(S) * V_T...
Subscript num_rows() const
BIAS::Vector< double > Upper_symmetric_eigenvalue_solve(const BIAS::Matrix< double > &A)
Solve symmetric eigenvalue problem (eigenvalues only)
void GetSystemMatrix(Matrix< T > &dest) const
compute square system matrix dest = A^T * A
const int RelRightNullspaceDim() const