Commit 61558c39 authored by Damien Leroux's avatar Damien Leroux
Browse files

Removed SVD from model::compute(). Too harmful. ColPiv* as well.

parent f140ec5a
......@@ -442,7 +442,7 @@ MatrixXd components(const MatrixXd& M, const MatrixXd& P)
}
enum class SolverType { QR, SVD };
enum class SolverType { QR };
typedef std::vector<std::map<model_block_key, MatrixXd>> constraint_list;
......@@ -575,22 +575,6 @@ struct model {
/*m_blocks.erase(m_blocks.begin() + (it - m_keys.begin()));*/
}
void use_SVD()
{
if (m_solver_type != SolverType::SVD) {
m_computed = false;
m_solver_type = SolverType::SVD;
}
}
void use_QR()
{
if (m_solver_type != SolverType::QR) {
m_computed = false;
m_solver_type = SolverType::QR;
}
}
void enable_constraints()
{
if (!m_with_constraints) {
......@@ -699,45 +683,28 @@ struct model {
* Either ColPivHouseholderQR or SVD give BAD results with our typical matrices.
* There is only FullPivHouseholderQR to give acceptable results.
*/
if (m_solver_type == SolverType::QR) {
/*MSG_DEBUG((*m_X));*/
int m_columns = m_X->outerSize();
/*ColPivHouseholderQR<MatrixXd> solver(*m_X);*/
FullPivHouseholderQR<MatrixXd> solver(*m_X);
{
MatrixXd Q = solver.matrixQ();
MatrixXd QR = solver.matrixQR();
/*MSG_DEBUG("Q" << std::endl << Q);*/
/*MSG_DEBUG("Q.Qt" << std::endl << (Q * Q.transpose()));*/
/*MSG_DEBUG("QR" << std::endl << QR);*/
}
/*MSG_DEBUG("Threshold: " << COMPONENT_EPSILON);*/
solver.setThreshold(COMPONENT_EPSILON);
m_rank = solver.rank();
/*MSG_DEBUG("Rank: " << m_rank);*/
m_coefficients->resize(m_columns, Y().outerSize());
/*MSG_DEBUG(MATRIX_SIZE((*m_coefficients)));*/
for (int i = 0; i < Y().outerSize(); ++i) {
/*MSG_DEBUG(MATRIX_SIZE(Y().col(i)) << std::endl << Y().col(i));*/
MatrixXd tmp = solver.solve(Y().col(i));
/*MSG_DEBUG(MATRIX_SIZE(tmp) << std::endl << tmp);*/
m_coefficients->col(i) = tmp;
}
} else {
MatrixXd x = *m_X;
JacobiSVD<MatrixXd> solver(x);
m_rank = 0;
int nzsv = solver.nonzeroSingularValues();
for (int i = 0; i < nzsv; ++i) {
m_rank += !around_zero(solver.singularValues()(i));
}
MatrixXd y = Y();
/* FIXME: there's a segfault somewhere around HERE. */
*m_coefficients = solver.solve(y);
}
/*MSG_DEBUG((*m_X));*/
int m_columns = m_X->outerSize();
FullPivHouseholderQR<MatrixXd> solver(*m_X);
{
MatrixXd Q = solver.matrixQ();
MatrixXd QR = solver.matrixQR();
/*MSG_DEBUG("Q" << std::endl << Q);*/
/*MSG_DEBUG("Q.Qt" << std::endl << (Q * Q.transpose()));*/
/*MSG_DEBUG("QR" << std::endl << QR);*/
}
/*MSG_DEBUG("Threshold: " << COMPONENT_EPSILON);*/
solver.setThreshold(COMPONENT_EPSILON);
m_rank = solver.rank();
/*MSG_DEBUG("Rank: " << m_rank);*/
m_coefficients->resize(m_columns, Y().outerSize());
/*MSG_DEBUG(MATRIX_SIZE((*m_coefficients)));*/
for (int i = 0; i < Y().outerSize(); ++i) {
/*MSG_DEBUG(MATRIX_SIZE(Y().col(i)) << std::endl << Y().col(i));*/
MatrixXd tmp = solver.solve(Y().col(i));
/*MSG_DEBUG(MATRIX_SIZE(tmp) << std::endl << tmp);*/
m_coefficients->col(i) = tmp;
}
/*MSG_DEBUG("Coefficients:" << std::endl << m_coefficients);*/
*m_residuals = Y() - X() * (*m_coefficients);
/*MSG_DEBUG("Residuals:" << std::endl << m_residuals);*/
......@@ -823,6 +790,7 @@ struct model {
MatrixXd XtX_pseudo_inverse() const
{
/* FIXME Should NOT use SVD (see note in compute() */
JacobiSVD<MatrixXd> inverter(XtX(), ComputeFullV);
auto& V = inverter.matrixV();
VectorXd inv_sv(inverter.singularValues());
......
......@@ -176,11 +176,6 @@ basic_model()
M.add_block({}, cross_indicator(all_populations()));
const char* decomp_method = getenv("DECOMP_METHOD");
std::string s(decomp_method ? decomp_method : "");
if (s == "SVD") {
M.use_SVD();
} else {
M.use_QR();
}
/*std::cout << "Xt" << std::endl << M.X().transpose() << std::endl;*/
/*std::cout << "Yt" << std::endl << M.Y().transpose() << std::endl;*/
/*std::cout << "XtX^-1" << std::endl << M.XtX_pseudo_inverse() << std::endl;*/
......@@ -381,11 +376,6 @@ init_model(const value<MatrixXd>& Y, const value<model_block_type>& indicator)
model M(Y, all_populations());
const char* decomp_method = getenv("DECOMP_METHOD");
std::string s(decomp_method ? decomp_method : "");
if (s == "SVD") {
M.use_SVD();
} else {
M.use_QR();
}
M.add_block({}, indicator);
/*std::cout << "indicator" << std::endl << indicator->transpose() << std::endl;*/
/*std::cout << "Xt" << std::endl << M.X().transpose() << std::endl;*/
......
Supports Markdown
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