|
|
Prev | Next |
| Matlab or Octave |
[y, r] = chol(x)
|
| C++ |
= cholesky(x, r)
|
= chol(x)
sets z to the
upper triangular @(@
m \times m
@)@ matrix such that
@[@
x = y^\T * y
@]@
If r is not zero, then the matrix x
is singular and r is the rank of the largest non-singular
principal minor.
function [ok] = cholesky_ok()
ok = true;
m = 3;
a = rand(m, m);
x = a' * a;
% -------------
[y, r] = chol(x);
% -------------
[m, n] = size(y);
ok = ok & (m == 3);
ok = ok & (n == 3);
ok = ok & all ( all( abs(y' * y - x) < 1e-10 ) );
for i = 1 : m
for j = 1 : (i-1)
ok = ok & (y(i,j) == 0);
end
end
ok = ok & (r == 0);
return
= cholesky(x, y, rank)
where x and y are
ublas matrix<double> objects.
The size of x and y are @(@
m \times m
@)@,
The return value of the size_t argument rank is the rank
of the largest principal minor of x.
# include <boost/numeric/ublas/matrix.hpp>
# include <mat2cpp.hpp>
bool cholesky_ok(void)
{ bool ok = true;
using namespace mat2cpp;
size_t i, j, m(3);
matrix<double> a = rand(m, m);
matrix<double> x = prod(trans(a), a);
// ---------------------------------------
size_t rank;
matrix<double> y = cholesky(x, rank);
// ---------------------------------------
ok &= (rank == m);
ok &= (y.size1() == m);
ok &= (y.size2() == m);
matrix<double> check = prod(trans(y), y);
for(i = 0; i < m; i++)
{ for(j = 0; j < m; j++)
ok &= std::fabs( check(i, j) - x(i,j) ) < 1e-10;
}
for(i = 0; i < m; i++)
{ for(j = 0; j < i; j++)
ok &= (y(i, j) == 0.);
}
return ok;
}