From 0f18eb836a66ee976d98e961b1954f4ced23f78d Mon Sep 17 00:00:00 2001
From: Mansard <nmansard@laas.fr>
Date: Tue, 23 Aug 2011 11:00:50 +0200
Subject: [PATCH] First complete working version.
---
unitTesting/qrt.cpp | 94 ++++++++++++++++++++++++++++++++++++---------
1 file changed, 75 insertions(+), 19 deletions(-)
diff --git a/unitTesting/qrt.cpp b/unitTesting/qrt.cpp
index 136864d..7792fec 100644
--- a/unitTesting/qrt.cpp
+++ b/unitTesting/qrt.cpp
@@ -23,36 +23,76 @@ std::ostream& p( const char * n ) { return std::cout << n << " = "; }
namespace Eigen
{
- class FullRowRankColPivQR
- : private ColPivHouseholderQR< MatrixXd >
+
+ /* This class add two methods to the classical ColPiv: solveInPlace for full
+ * col rank matrices (and the direct consequence of pseudoInverse for
+ * f-c-r). And solveTranposeInPlace for full row rank matrices (and the
+ * direct pseudoInverseTranspose).
+ */
+ class ColPivQRSolveInPlace
+ : public ColPivHouseholderQR< MatrixXd >
{
public:
- void compute(const MatrixType& matrix)
+ /* Find a solution x to the problem G x = b. The input vector is b,
+ * completed by the 0 tail to obtain the x size. The solver only works for
+ * full-col rank matrices, ie matrices G = Q [ R; 0 ], with R full-rank
+ * (full diag) upper triangular. */
+ void solveInPlace( MatrixXd& Gp )
{
- assert( matrix.rows()<=matrix.cols() );
- ColPivHouseholderQR< MatrixXd >::compute( matrix.transpose() );
- assert( rank() == matrix.rows() );
+ const int r = rank(), m=cols(), n=rows();
+ assert( r==m );
+ assert( Gp.rows() == m ); // TODO: if not proper size, resize.
+
+ VectorXd workspace( n );
+ /* P2*P1*P0 ... */
+ for (int k = r-1; k >= 0 ; --k)
+ {
+ int remainingSize = n-k;
+ Gp.rightCols( Gp.cols()-k )
+ .applyHouseholderOnTheRight(matrixQR().col(k).tail(remainingSize-1),
+ hCoeffs().coeff(k), &workspace.coeffRef(0));
+ }
+
+ matrixQR()
+ .topLeftCorner(r,r).triangularView<Upper>()
+ .solveInPlace(Gp);
+
+ Gp.applyOnTheLeft( colsPermutation() );
}
- void solveInPlace( MatrixXd& Gp )
+ MatrixXd pseudoInverse( void )
{
- /* r is the rank, nxm the size of the original matrix (ie whose transpose has
- * been decomposed. */
+ MatrixXd res = MatrixXd::Identity( cols(),rows() );
+ solveInPlace( res );
+ return res;
+ }
+
+
+
+ /* Find a solution x to the problem G'x = b. The input vector is b,
+ * completed by the 0 tail to obtain the x size. The solver only works for
+ * full-row rank matrices, ie matrices G' = [ L 0 ] Q, with L full-rank
+ * (full diag) lower triangular. */
+ void solveTransposeInPlace( MatrixXd& Gtp )
+ {
+ /* r is the rank, nxm the size of the original matrix (ie whose transpose
+ * has been decomposed. n is the number of cols of the original matrix,
+ * thus number of rows of the transpose we want to inverse. */
const int r = rank(), n=cols(), m=rows();
assert( r==n );
- assert( Gp.rows() == m );
+ assert( Gtp.rows() == m ); // TODO: if not proper size, resize.
- /* G'E = Q R
- * => G^+ = Q R^+T E' */
+ /* G E = Q R :: E' G' = L Q', with L=R'
+ * => G^+T = Q R^+T E' */
/* Compute E'*X. */
- Gp.applyOnTheLeft( colsPermutation().transpose() );
+ Gtp.applyOnTheLeft( colsPermutation().transpose() );
/* Compute R^+T E'*X. */
matrixQR()
.topLeftCorner(r,r).transpose().triangularView<Lower>()
- .solveInPlace(Gp.topRows(r));
+ .solveInPlace(Gtp.topRows(r));
/* Compute Q R^+T E'*X. */
/* Q = P1*P2* ... *Pn */
@@ -60,11 +100,19 @@ namespace Eigen
for (int k = r-1; k >= 0 ; --k)
{
int remainingSize = m-k;
- Gp.bottomRows( Gp.rows()-k )
+ Gtp.bottomRows( Gtp.rows()-k )
.applyHouseholderOnTheLeft(matrixQR().col(k).tail(remainingSize-1),
hCoeffs().coeff(k), &workspace.coeffRef(0));
}
}
+
+ MatrixXd pseudoInverseTranspose( void )
+ {
+ MatrixXd Gtp = MatrixXd::Identity( rows(),cols() );
+ solveTransposeInPlace( Gtp );
+ return Gtp;
+ }
+
};
}
@@ -75,7 +123,7 @@ int main (void)
{ /* Compute the pseudo inverse of a full-col rank matrix: G = V [ R; 0 ]. */
- const unsigned int n=6, m=3, r=3;
+ const unsigned int n=12, m=9, r=9;
MatrixXd G = MatrixXd::Random(n,m);
p("G") << (MATLAB)G << endl;
@@ -134,10 +182,17 @@ int main (void)
Gp.applyOnTheLeft( G_qr.colsPermutation() );
p("Gp") << (MATLAB)Gp << endl;
}
+
+ {
+ ColPivQRSolveInPlace G_qr; G_qr.compute(G);
+ p("Gp") << (MATLAB)G_qr.pseudoInverse() << endl;
+ }
+
+
}
{ /* Compute the pseudo inverse of a full-row rank matrix: G = [ L; 0 ] V'. */
- const unsigned int n=3, m=6, r=3;
+ const unsigned int n=9, m=12, r=9;
MatrixXd G = MatrixXd::Random(n,m);
p("G") << (MATLAB)G << endl;
@@ -174,10 +229,11 @@ int main (void)
}
{
- FullRowRankColPivQR G_lq; G_lq.compute(G);
+ ColPivQRSolveInPlace Gt_qr; Gt_qr.compute(G.transpose());
MatrixXd Gp = MatrixXd::Identity( m,r );
- G_lq.solveInPlace(Gp);
+ Gt_qr.solveTransposeInPlace(Gp);
p("Gp") << (MATLAB)Gp << endl;
+ p("Gp") << (MATLAB)Gt_qr.pseudoInverseTranspose() << endl;
}
}
--
GitLab