Commit 00197638 authored by Davis King's avatar Davis King
Browse files

Changed the pinv() function so that it doesn't get slow when operating on

matrices with many more columns than rows.

--HG--
extra : convert_revision : svn%3Afdd8eb12-d10e-0410-9acb-85c331704f74/trunk%403029
parent d65a1009
...@@ -1278,9 +1278,15 @@ convergence: ...@@ -1278,9 +1278,15 @@ convergence:
template < template <
typename EXP typename EXP
> >
inline const matrix<typename EXP::type,EXP::NC,EXP::NR,typename EXP::mem_manager_type> pinv ( const matrix<typename EXP::type,EXP::NC,EXP::NR,typename EXP::mem_manager_type> pinv_helper (
const matrix_exp<EXP>& m const matrix_exp<EXP>& m
) )
/*!
ensures
- computes the results of pinv(m) but does so using a method that is fastest
when m.nc() <= m.nr(). So if m.nc() > m.nr() then it is best to use
trans(pinv_helper(trans(m))) to compute pinv(m).
!*/
{ {
typename matrix_exp<EXP>::matrix_type u; typename matrix_exp<EXP>::matrix_type u;
typedef typename EXP::mem_manager_type MM1; typedef typename EXP::mem_manager_type MM1;
...@@ -1307,6 +1313,21 @@ convergence: ...@@ -1307,6 +1313,21 @@ convergence:
return tmp(scale_columns(v,reciprocal(round_zeros(w,eps))))*trans(u); return tmp(scale_columns(v,reciprocal(round_zeros(w,eps))))*trans(u);
} }
template <
typename EXP
>
const matrix<typename EXP::type,EXP::NC,EXP::NR,typename EXP::mem_manager_type> pinv (
const matrix_exp<EXP>& m
)
{
// if m has more columns then rows then it is more efficient to
// compute the pseudo-inverse of its transpose (given the way I'm doing it below).
if (m.nc() > m.nr())
return trans(pinv_helper(trans(m)));
else
return pinv_helper(m);
}
// ---------------------------------------------------------------------------------------- // ----------------------------------------------------------------------------------------
template < template <
......
...@@ -121,7 +121,6 @@ namespace ...@@ -121,7 +121,6 @@ namespace
DLIB_TEST((equal(round_zeros(mi*m,0.000001) , identity_matrix<double,5>()))); DLIB_TEST((equal(round_zeros(mi*m,0.000001) , identity_matrix<double,5>())));
} }
{ {
matrix<double,5,2,MM,column_major_layout> m; matrix<double,5,2,MM,column_major_layout> m;
...@@ -162,6 +161,46 @@ namespace ...@@ -162,6 +161,46 @@ namespace
DLIB_TEST((equal(round_zeros(mi*m,0.000001) , identity_matrix<double,2>()))); DLIB_TEST((equal(round_zeros(mi*m,0.000001) , identity_matrix<double,2>())));
} }
{
matrix<double,5,2,MM,column_major_layout> m;
for (long r = 0; r < m.nr(); ++r)
{
for (long c = 0; c < m.nc(); ++c)
{
m(r,c) = r*c;
}
}
m = cos(exp(m));
matrix<double> mi = trans(pinv(trans(m) ));
DLIB_TEST(mi.nr() == m.nc());
DLIB_TEST(mi.nc() == m.nr());
DLIB_TEST((equal(round_zeros(mi*m,0.000001) , identity_matrix<double,2>())));
}
{
matrix<double> m(5,2);
for (long r = 0; r < m.nr(); ++r)
{
for (long c = 0; c < m.nc(); ++c)
{
m(r,c) = r*c;
}
}
m = cos(exp(m));
matrix<double> mi = trans(pinv(trans(m) ));
DLIB_TEST(mi.nr() == m.nc());
DLIB_TEST(mi.nc() == m.nr());
DLIB_TEST((equal(round_zeros(mi*m,0.000001) , identity_matrix<double,2>())));
}
{ {
matrix<long> a1(5,1); matrix<long> a1(5,1);
matrix<long,0,0,MM,column_major_layout> a2(1,5); matrix<long,0,0,MM,column_major_layout> a2(1,5);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment