refs#16

In NumA++ (trunk):
- New LSMR ported from scipy.

include/NumA/linear_algebra/LinAlgUtils.h

@@ -8,9 +8,15 @@
 #ifndef LINALGUTILS_H_
 #define LINALGUTILS_H_
 
+#include <cstddef>
+
 #include "matrix/MatrixD.h"
 #include "vector/VectorD.h"
 
+namespace ElemUtils {
+class Formatter;
+} /* namespace ElemUtils */
+
 namespace NumA {
 class MatrixD;
 } /* namespace NumA */
@@ -38,6 +44,17 @@ public:
     static VectorD solve(const MatrixD & A, const VectorD & Y, Method method =
             ColQR);
 
+    /**
+     * @brief Solves the linear system \f$ A X = Y \f$ and returns \f$ X \f$ in the least-squares sense.
+     * Iterative method for large sparse matrices.
+     * @param A MatrixD
+     * @param Y VectorD
+     * @return VectorD Solution X.
+     */
+    static VectorD solveLSMR(const MatrixD & A, const VectorD & Y, double damp =
+            0., double atol = 1.e-6, double btol = 1.e-6, double conlim = 1.e8,
+            size_t maxiter = 0, ElemUtils::Formatter* output = 0);
+
     /**
      * @brief Computes the rank of the matrix with the given method.
      * @param A MatrixD
@@ -45,6 +62,9 @@ public:
      * @return rank unsigned int
      */
     static unsigned int rank(const MatrixD & A, Method method = ColQR);
+
+private:
+    static void _sym_ortho(double a, double b, double& c, double& s, double& r);
 };
 
 } // namespace NumA

include/NumA/linear_algebra/vector/VectorD.h

@@ -104,17 +104,29 @@ public:
 
     /**
      * Scalar product
-     * @param rhs
+     * @param rhs VectorD of same size
      * @return double
      */
     double operator *(const VectorD &rhs) const;
 
     /**
      * Subtraction.
-     * @param rhs
+     * @param rhs VectorD of same size
      * @return VectorD of same size
      */
     VectorD operator -(const VectorD &rhs) const;
+    /**
+     * Addition
+     * @param rhs VectorD of same size
+     * @return VectorD of same size
+     */
+    VectorD operator +(const VectorD &rhs) const;
+
+    // Vector/double operations
+    VectorD operator*(double rhs) const;
+    VectorD operator+(double rhs) const;
+    VectorD operator-(double rhs) const;
+    VectorD operator/(double rhs) const;
 
     /**
      * Norm 2.

src/NumA/linear_algebra/LinAlgUtils.cpp

 #include "../../../include/NumA/linear_algebra/LinAlgUtils.h"
 
+#include <algorithm>
+#include <cmath>
+#include <ctime>
 #include <Eigen/Core>
 #include <Eigen/LU>
 #include <Eigen/QR>
-#include <sstream>
-#include <stdexcept>
+#include <limits>
+#include <string>
+#include <vector>
+#include <ElementaryUtils/logger/CustomException.h>
+#include <ElementaryUtils/string_utils/Formatter.h>
 
 #include "../../../include/NumA/linear_algebra/eigen/EigenUtils.h"
 
@@ -13,16 +19,18 @@ namespace NumA {
 VectorD LinAlgUtils::solve(const MatrixD & A, const VectorD & Y,
         Method method) {
     if (A.cols() != A.rows()) {
-        std::stringstream formatter;
-        formatter << "[Solve::solve] Matrix column size = " << A.cols()
+        ElemUtils::Formatter formatter;
+        formatter << "Matrix column size = " << A.cols()
                 << " does not match matrix row size = " << A.rows();
-        throw std::runtime_error(formatter.str());
+        throw ElemUtils::CustomException("LinAlgUtils", __func__,
+                formatter.str());
     }
     if (Y.size() != A.rows()) {
-        std::stringstream formatter;
-        formatter << "[Solve::solve] Vector's size = " << Y.size()
+        ElemUtils::Formatter formatter;
+        formatter << "Vector's size = " << Y.size()
                 << " does not match matrix size = " << A.rows();
-        throw std::runtime_error(formatter.str());
+        throw ElemUtils::CustomException("LinAlgUtils", __func__,
+                formatter.str());
     }
 
     // Converting matrix and vector to Eigen
@@ -55,6 +63,339 @@ VectorD LinAlgUtils::solve(const MatrixD & A, const VectorD & Y,
     return X;
 }
 
+VectorD LinAlgUtils::solveLSMR(const MatrixD & A, const VectorD & Y,
+        double damp, double atol, double btol, double conlim, size_t maxiter,
+        ElemUtils::Formatter* output) {
+
+    std::vector<std::string> msg(8);
+    msg[0] = "The exact solution is x = 0 ";
+    msg[1] = "Ax - b is small enough, given atol, btol ";
+    msg[2] = "The least-squares solution is good enough, given atol ";
+    msg[3] = "The estimate of cond(Abar) has exceeded conlim ";
+    msg[4] = "Ax - b is small enough for this machine ";
+    msg[5] = "The least-squares solution is good enough for this machine";
+    msg[6] = "Cond(Abar) seems to be too large for this machine ";
+    msg[7] = "The iteration limit has been reached ";
+
+    std::string hdg1 = "  itn x(1) norm r norm A'r";
+    std::string hdg2 = " compatible LS norm A cond A";
+    unsigned int pfreq = 20;   // print frequency (for repeating the heading)
+    size_t pcount = 0;  // print counter
+
+    // Converting matrix and vector to Eigen
+    Eigen::MatrixXd EigenA = EigenUtils::convertToEigen(A);
+    Eigen::VectorXd EigenY = EigenUtils::convertToEigen(Y);
+
+    size_t m = A.rows();
+    size_t n = A.cols();
+//    MatrixD At = A.transpose();
+    Eigen::MatrixXd EigenAt = EigenA.transpose();
+
+    // stores the num of singular values
+    size_t minDim = std::min(m, n);
+
+    if (maxiter <= 0)
+        maxiter = minDim;
+
+    if (output) {
+        *output << " ";
+        *output << "LSMR Least-squares solution of Ax = b\n";
+        *output << "The matrix A has " << m << " rows and " << n << " cols\n";
+        *output << "damp = " << damp << "\n";
+        *output << "atol = " << atol << "  conlim = " << conlim << "\n";
+        *output << "btol = " << btol << "  maxiter = " << maxiter << "\n";
+    }
+
+//    VectorD u(Y);
+    Eigen::VectorXd u(EigenY);
+    double beta = u.norm();
+
+//    VectorD v(n);
+    Eigen::VectorXd v(n);
+    double alpha = 0.;
+
+    if (beta > 0) {
+        u = u / beta;
+//        v = At * u;
+        v = EigenAt * u;
+        alpha = v.norm();
+    }
+
+    if (alpha > 0) {
+        v = v / alpha;
+    }
+
+    // Initialize variables for 1st iteration.
+
+    size_t itn = 0;
+    double zetabar = alpha * beta;
+    double alphabar = alpha;
+    double rho = 1;
+    double rhobar = 1;
+    double cbar = 1;
+    double sbar = 0;
+
+//    VectorD h(v);
+//    VectorD hbar(n);
+//    VectorD x(n);
+    Eigen::VectorXd h(v);
+    Eigen::VectorXd hbar(n);
+    Eigen::VectorXd x(n);
+
+    // Initialize variables for estimation of ||r||.
+
+    double betadd = beta;
+    double betad = 0.;
+    double rhodold = 1.;
+    double tautildeold = 0.;
+    double thetatilde = 0.;
+    double zeta = 0.;
+    double d = 0.;
+
+    // Initialize variables for estimation of ||A|| and cond(A)
+
+    double normA2 = alpha * alpha;
+    double maxrbar = 0.;
+    double minrbar = 1.e+100;
+    double normA = std::sqrt(normA2);
+    double condA = 1.;
+    double normx = 0.;
+
+    // Items for use in stopping rules.
+    double normb = beta;
+    size_t istop = 0;
+    double ctol = 0.;
+    if (conlim > 0.)
+        ctol = 1. / conlim;
+    double normr = beta;
+
+    // Reverse the order here from the original matlab code because
+    // there was an error on return when arnorm==0
+    double normar = alpha * beta;
+    if (normar == 0.) {
+        if (output)
+            *output << msg[0];
+        return EigenUtils::convertToVectorD(x); //, istop, itn, normr, normar, normA, condA, normx
+    }
+
+    double test1, test2, test3, t1, rtol;
+    if (output) {
+        *output << "\n";
+        *output << hdg1 << " " << hdg2 << "\n";
+        test1 = 1.;
+        test2 = alpha / beta;
+        *output << itn << " " << x[0] << " " << normr << " " << normar << " "
+                << test1 << " " << test2 << "\n";
+    }
+
+    // Initialize variables used in the loop
+    double chat, shat, alphahat, rhoold, c, s, thetanew, rhobarold, zetaold,
+            thetabar, rhotemp, betaacute, betacheck, betahat, thetatildeold,
+            ctildeold, stildeold, rhotildeold, taud, betadtaud;
+
+    // Test time of an iteration
+    std::clock_t start;
+    double duration;
+
+    // Main iteration loop.
+    while (itn < maxiter) {
+        itn++;
+
+        // Test time of an iteration
+        start = std::clock();
+
+        // Perform the next step of the bidiagonalization to obtain the
+        // next  beta, u, alpha, v.  These satisfy the relations
+        //         beta*u  =  a*v   -  alpha*u,
+        //        alpha*v  =  A'*u  -  beta*v.
+
+        u = (EigenA * v) - (u * alpha);
+        beta = u.norm();
+
+        if (beta > 0) {
+            u = u / beta;
+//            v = (At * u) - (v * beta);
+            v = (EigenAt * u) - (v * beta);
+            alpha = v.norm();
+            if (alpha > 0)
+                v = v / alpha;
+        }
+
+        // At this point, beta = beta_{k+1}, alpha = alpha_{k+1}.
+
+        // Construct rotation Qhat_{k,2k+1}.
+
+        _sym_ortho(alphabar, damp, chat, shat, alphahat);
+
+        // Use a plane rotation (Q_i) to turn B_i to R_i
+
+        rhoold = rho;
+        _sym_ortho(alphahat, beta, c, s, rho);
+        thetanew = s * alpha;
+        alphabar = c * alpha;
+
+        // Use a plane rotation (Qbar_i) to turn R_i^T to R_i^bar
+
+        rhobarold = rhobar;
+        zetaold = zeta;
+        thetabar = sbar * rho;
+        rhotemp = cbar * rho;
+        _sym_ortho(cbar * rho, thetanew, cbar, sbar, rhobar);
+        zeta = cbar * zetabar;
+        zetabar = -sbar * zetabar;
+
+        // Update h, h_hat, x.
+
+        hbar = h - (hbar * (thetabar * rho / (rhoold * rhobarold)));
+        x = x + (hbar * (zeta / (rho * rhobar)));
+        h = v - (h * (thetanew / rho));
+
+        // Estimate of ||r||.
+
+        // Apply rotation Qhat_{k,2k+1}.
+        betaacute = chat * betadd;
+        betacheck = -shat * betadd;
+
+        // Apply rotation Q_{k,k+1}.
+        betahat = c * betaacute;
+        betadd = -s * betaacute;
+
+        // Apply rotation Qtilde_{k-1}.
+        // betad = betad_{k-1} here.
+
+        thetatildeold = thetatilde;
+        _sym_ortho(rhodold, thetabar, ctildeold, stildeold, rhotildeold);
+        thetatilde = stildeold * rhobar;
+        rhodold = ctildeold * rhobar;
+        betad = -stildeold * betad + ctildeold * betahat;
+
+        // betad   = betad_k here.
+        // rhodold = rhod_k  here.
+
+        tautildeold = (zetaold - thetatildeold * tautildeold) / rhotildeold;
+        taud = (zeta - thetatilde * tautildeold) / rhodold;
+        d = d + betacheck * betacheck;
+        betadtaud = betad - taud;
+        normr = std::sqrt(d + betadtaud * betadtaud + betadd * betadd);
+
+        // Estimate ||A||.
+        normA2 = normA2 + beta * beta;
+        normA = std::sqrt(normA2);
+        normA2 = normA2 + alpha * alpha;
+
+        // Estimate cond(A).
+        maxrbar = std::max(maxrbar, rhobarold);
+        if (itn > 1)
+            minrbar = std::min(minrbar, rhobarold);
+        condA = std::max(maxrbar, rhotemp) / std::min(minrbar, rhotemp);
+
+        // Test for convergence.
+
+        // Compute norms for convergence testing.
+        normar = std::abs(zetabar);
+        normx = x.norm();
+
+        // Now use these norms to estimate certain other quantities,
+        // some of which will be small near a solution.
+
+        test1 = normr / normb;
+        if ((normA * normr) != 0.)
+            test2 = normar / (normA * normr);
+        else
+            test2 = std::numeric_limits<double>::infinity();
+        test3 = 1 / condA;
+        t1 = test1 / (1 + normA * normx / normb);
+        rtol = btol + atol * normA * normx / normb;
+
+        // The following tests guard against extremely small values of
+        // atol, btol or ctol.  (The user may have set any or all of
+        // the parameters atol, btol, conlim  to 0.)
+        // The effect is equivalent to the normAl tests using
+        // atol = eps,  btol = eps,  conlim = 1/eps.
+
+        if (itn >= maxiter)
+            istop = 7;
+        if (1 + test3 <= 1)
+            istop = 6;
+        if (1 + test2 <= 1)
+            istop = 5;
+        if (1 + t1 <= 1)
+            istop = 4;
+
+        // Allow for tolerances set by the user.
+
+        if (test3 <= ctol)
+            istop = 3;
+        if (test2 <= atol)
+            istop = 2;
+        if (test1 <= rtol)
+            istop = 1;
+
+        // Test time of an iteration
+        duration = (std::clock() - start) / (double) CLOCKS_PER_SEC;
+
+        // See if it is time to print something.
+
+        if (output) {
+            if ((n <= 40) or (itn <= 10) or (itn >= maxiter - 10)
+                    or (itn % 10 == 0) or (test3 <= 1.1 * ctol)
+                    or (test2 <= 1.1 * atol) or (test1 <= 1.1 * rtol)
+                    or (istop != 0)) {
+                if (pcount >= pfreq) {
+                    pcount = 0;
+                    *output << "\n";
+                    *output << hdg1 << hdg2 << "\n";
+                }
+                pcount = pcount + 1;
+                *output << itn << " " << x[0] << " " << normr << " " << normar
+                        << " " << test1 << " " << test2 << " " << normA << " "
+                        << condA << " " << duration << "s" << "\n";
+            }
+        }
+
+        if (istop > 0)
+            break;
+
+    }
+
+    // Print the stopping condition.
+
+    if (output) {
+        *output << "\n";
+        *output << "LSMR finished" << "\n";
+        *output << msg[istop] << "\n";
+        *output << "istop = " << istop << "  normr = " << normr << "\n";
+        *output << "  normA = " << normA << " normAr = " << normar << "\n";
+        *output << "itn = " << itn << " condA = " << condA << "\n";
+        *output << "  normx =" << normx << "\n";
+    }
+
+    return EigenUtils::convertToVectorD(x); //, istop, itn, normr, normar, normA, condA, normx
+}
+
+void LinAlgUtils::_sym_ortho(double a, double b, double& c, double& s,
+        double& r) {
+    if (b == 0) {
+        c = (a > 0) - (a < 0);
+        s = 0;
+        r = std::abs(a);
+    } else if (a == 0) {
+        c = 0;
+        s = (b > 0) - (b < 0);
+        r = std::abs(b);
+    } else if (std::abs(b) > std::abs(a)) {
+        double tau = a / b;
+        s = ((b > 0) - (b < 0)) / std::sqrt(1. + tau * tau);
+        c = s * tau;
+        r = b / s;
+    } else {
+        double tau = b / a;
+        c = ((a > 0) - (a < 0)) / std::sqrt(1. + tau * tau);
+        s = c * tau;
+        r = a / c;
+    }
+}
+
 unsigned int LinAlgUtils::rank(const MatrixD& A, Method method) {
     // Converting matrix and vector to Eigen
     Eigen::MatrixXd EigenA = EigenUtils::convertToEigen(A);

src/NumA/linear_algebra/matrix/MatrixD.cpp

@@ -348,7 +348,7 @@ const double& MatrixD::at(const size_t i, const size_t j) const {
                 << rows();
         throw std::runtime_error(formatter.str());
     }
-    return m_matrix.at(i * cols() + j);
+    return m_matrix[i * cols() + j];
 }
 
 } /* namespace NumA */

src/NumA/linear_algebra/vector/VectorD.cpp

@@ -2,10 +2,11 @@
 
 #include <cmath>
 #include <sstream>
-#include <stdexcept>
+//#include <stdexcept>
+#include <ElementaryUtils/logger/CustomException.h>
 
-#include "../../../../include/NumA/linear_algebra/vector/Vector4D.h"
 #include "../../../../include/NumA/linear_algebra/matrix/MatrixD.h"
+#include "../../../../include/NumA/linear_algebra/vector/Vector4D.h"
 
 namespace NumA {
 
@@ -77,7 +78,7 @@ double VectorD::operator *(const VectorD &rhs) const {
                 "[VectorD::operator *] Vectors have different size");
     }
 
-    for (unsigned int i = 0; i != vectorSize; i++) {
+    for (unsigned int i = 0; i < vectorSize; i++) {
         result += (m_vector[i] * rhs[i]);
     }
 
@@ -92,12 +93,58 @@ VectorD NumA::VectorD::operator -(const VectorD& rhs) const {
     }
 
     VectorD result(size());
-    for (unsigned int i = 0; i != size(); i++) {
+    for (unsigned int i = 0; i < size(); i++) {
         result[i] = (m_vector[i] - rhs[i]);
     }
     return result;
 }
 
+VectorD VectorD::operator +(const VectorD& rhs) const {
+    // cannot perform addition if the two vector's size are different
+    if (size() != rhs.size()) {
+        throw ElemUtils::CustomException("VectorD", __func__,
+                "Vectors have different size.");
+    }
+
+    VectorD result(size());
+    for (unsigned int i = 0; i < size(); i++) {
+        result[i] = (m_vector[i] + rhs[i]);
+    }
+    return result;
+}
+
+VectorD VectorD::operator *(double rhs) const {
+    VectorD result(size());
+    for (unsigned int i = 0; i < size(); i++) {
+        result[i] = (m_vector[i] * rhs);
+    }
+    return result;
+}
+
+VectorD VectorD::operator +(double rhs) const {
+    VectorD result(size());
+    for (unsigned int i = 0; i < size(); i++) {
+        result[i] = (m_vector[i] + rhs);
+    }
+    return result;
+}
+
+VectorD VectorD::operator -(double rhs) const {
+    VectorD result(size());
+    for (unsigned int i = 0; i < size(); i++) {
+        result[i] = (m_vector[i] - rhs);
+    }
+    return result;
+}
+
+VectorD VectorD::operator /(double rhs) const {
+    VectorD result(size());
+    for (unsigned int i = 0; i < size(); i++) {
+        result[i] = (m_vector[i] / rhs);
+    }
+    return result;
+}
+
 double NumA::VectorD::norm() const {
     double result = (*this) * (*this);
     result = std::sqrt(result);
@@ -124,7 +171,7 @@ VectorD VectorD::sub(size_t startPos, size_t endPos) const {
     VectorD result(static_cast<int>(endPos) - static_cast<int>(startPos));
 
     for (size_t i = 0; i < result.size(); i++) {
-        result[i]= m_vector[i+startPos];
+        result[i] = m_vector[i + startPos];
     }
 
     return result;