[eigen] OneNorm and condition number estimation
• To: eigen@xxxxxxxxxxxxxxxxxxx
• Subject: [eigen] OneNorm and condition number estimation
• From: R Kannan <ltbuckling+eigen@xxxxxxxxx>
• Date: Fri, 13 Jan 2012 18:22:58 +0000
• Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=mime-version:sender:date:x-google-sender-auth:message-id:subject :from:to:content-type; bh=SOFE6/Zu94U+BDErs2UAHCEsDshKnA2TF6G75GAbz5M=; b=IdkjrW9JvOgKtJrVCQY2myuSfGTBcTVeZkHesBahBj6GSJW49EDwb1/eHQKJWqYGLk DClaydPBr+0CuyViiGPEd95dmu6wbfwIAs8Fdqrvco3P48fh0YyVwSyxGfK2JCQyvGY/ 2W7vQirZmu0apF9wX+0Z4gZ8J0XImgHj3MoWk=

All,

I've just finished implementing an Eigen one norm estimator based on Higham and Tisseur's algorithm[1]. (attached code below.)

The OneNormEst class works similar to the Matlab command normest1. It doesn't need to know about A, it only needs 2 functors to evaluate Ax and A^T x.

EigenSparseMatrix sm(aux); // symmetric

Matrix::OneNormEst estimator(sm.rows(), 4);
double norm_a;

auto applyA = [&](EigenDenseMatrix& B, EigenDenseMatrix& Y)
{
};

estimator.ANorm(applyA, applyA, norm_a);

In conjunction with a functor that can solve Ax=b, it can be used to evaluate a condition number estimate, where the condition number kappa(A) is defined as
kappa(A) = ||A||_1 * ||A^{-1}||_1.
(which is what Matlab's condest does).

I'm happy to submit a patch -- does anyone have any views on where it can go? The code makes liberal use of lambdas so it's unsuitable for C++98.

[1] Nicholas J. Higham and Fran\c{c}oise Tisseur, A Block Algorithm for Matrix 1-Norm Estimation with an  Application to 1-Norm Pseudospectra, SIAM J. Matrix Anal. App. 21, 1185-1201, 2000.

-------code--------
#include <Eigen/Core>
#include <Eigen/Sparse>

#include <cmath>
#include <random>
#include <vector>

using std::cout;
static std::mt19937 engine;
static const int itmax=5;

namespace Matrix

{

class OneNormEst
{
public:
OneNormEst(const int& n, const int& t=4): m_n(n), m_t(t) {}
~OneNormEst() {}

template<typename A_operator, typename A_transpose_operator>
bool ANorm(const A_operator& applyA, A_transpose_operator& applyA_trans, double& norm)
{
Eigen::MatrixXd X(m_n, m_t), Y(m_n, m_t), Z(m_n, m_t);

std::uniform_int_distribution<int> rand_gen(0, m_n-1);
auto random_plus_minus_1_func = [&](const double& ) -> double
{
if (rand_gen(engine)>m_n/2)
return 1;
else
return -1;
};

X = X.unaryExpr(random_plus_minus_1_func);
X.col(0).setOnes();
X /= m_n;

Eigen::MatrixXd S(m_n, m_t), S_old(m_n, m_t); // sign matrix
Eigen::MatrixXi prodS(m_t, m_t);
S.setZero();

double est=0., est_old = 0.;

int ind_best;
std::vector<int> indices(m_n);
std::vector<int> ind_hist;

for (int k(0); k<itmax; ++k)
{
applyA(X, Y); // Y = A * X

int ind(-1);
for(int i(0); i<m_t; ++i)
{
double norm = Y.col(i).cwiseAbs().sum(); // norm = {||Y||}_1
if (norm > est) {
est = norm;
ind = indices[i];
}
}

if (est > est_old || k==1)
ind_best = ind;

if (est < est_old && k >= 1)
{
norm = est_old;
return true;
}

est_old = est;
S_old = S;

// S = sign(Y)
S = Y.unaryExpr([&] (const double& coeff) -> double
{
if(coeff >= 0.) return 1.;
else            return -1.;
});

prodS = (S_old.transpose() * S).matrix().cast<int>();

// if all cols are parallel, prodS will have all entries = n
if (prodS.cwiseAbs().colwise().sum().sum() == m_n*m_n*m_t) {
norm = est;
return true;
}

// if S_old(i) is parallel to S(j), replace S(j) with random
for(int i(0); i<m_t; ++i)
{
for(int j(0); j<m_t; ++j)
{
if(prodS.coeff(i, j)==m_n)
S.col(j) = S.col(j).unaryExpr(random_plus_minus_1_func);
}
}

// if S(i) is parallel to S(j), replace S(j) with random
prodS = (S.transpose()*S).matrix().cast<int>();

for(int i(0); i<m_t; ++i)
{
for(int j(i+1); j<m_t; ++j)
{
if(prodS.coeff(i, j)==m_n)
S.col(j) = S.col(j).unaryExpr(random_plus_minus_1_func);
}
}

applyA_trans(S, Z);//Z = A_transpose * S

Eigen::VectorXd h = Z.cwiseAbs().rowwise().maxCoeff();

if (k >= 1 && h.maxCoeff() == h.coeff(ind_best) ) {
norm = est;
return true;
}

for(int i(0); i<m_n; ++i)
indices[i] = i;

// sort indices by entries in h
std::sort(indices.begin(), indices.end(), [&](int& left, int& right)
{
return h.coeff(left) > h.coeff(right);
});

int n_found(0);
for(auto i=indices.begin(); i!=indices.begin()+m_t; ++i) // is indices contained in ind_hist?
if(std::find(ind_hist.begin(), ind_hist.end(), *i) != ind_hist.end())
++n_found;

if(n_found == m_t) {
norm = est;
return true;
}

std::vector<int> new_indices;

if(k>0)
{
new_indices.reserve(m_t);
int count(0);
for(auto it=indices.begin()+m_t; it !=indices.end() && count <m_t; ++it)
{
if(std::find(ind_hist.begin(), ind_hist.end(), *it) == ind_hist.end()) {
new_indices.push_back(*it); ++count;
}
}
new_indices.swap(indices);
}
assert(indices.size()>0);

X.setZero();
for(auto i=0; i<m_t; ++i)
{
X.coeffRef(indices[i], i) = 1; //Set X(:, j) = e_ind-j

ind_hist.push_back(indices[i]);
}
}
norm = est;

return true;

}

private:
int m_n; // rows
int m_t; // cols

};

};

 Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/