Hi, I'm trying to implement some form of Kaczmarz iterations. My code is as follows:

Code: Select all

typedef Eigen::Matrix<long double, Eigen::Dynamic, 1> VectorXld;
VectorXld Kaczmarz(const Eigen::SparseMatrix<long double> &A, const VectorXld &b, const VectorXld &x0) {
      long double d = 1;
      int row = 0, steps = 0;
      VectorXld x = x0;
      Eigen::SparseMatrix<long double> rowA;
      Eigen::SparseMatrix<long double> e;
      e.resize(1, A.rows());
      VectorXld rem = A * x - b;
      long double b_norm = b.norm();
      while ((d > 1e-5) || (steps < A.rows())) {

         // constructs a row vector (0, 0, ..., 0, 1, 0, ..., 0) with 1 at position "row"
         e.data().squeeze();
         e.reserve(1);
         e.coeffRef(0, row) = 1;

         // extracts "row"-th row from the matrix A
         rowA = e * A;

         // ( b_i - < a_i, x^(k) > ) / || a_i ||^2
         long double factor = (b.coeff(row, 0) - (rowA*x).coeff(0, 0)) / rowA.squaredNorm();

         // x^(k+1) = x^(k) + (...) * a_i^(T)
         for (int i = 0; i < x.rows(); i++) {
            x.coeffRef(i, 0) += factor * rowA.coeff(0, i);
         }
         rem = A * x - b;
         d = rem.norm() / b_norm;
         row++; if (row == A.rows()) row = 0;
         steps++;
      }
      return x;
   }

The problem is: this is quite slow. I think that the bottleneck can be nailed down by answering these questions:

1. is there a more (ideally a WAY more) efficient to extract a sparse row from a sparse matrix, other than e * A (as I'm doing it in my code from the lack of any other idea)?
2. is there a more suitable way to get the dot product between two vectors than (x*y).coeff(0, 0)?
3. what is the cost of transposition of a row vector into a column vector? What is the cost of transposition of a sparse matrix?
4. how would you effectively compute the expression: s^T A A^T s, where s is a column vector (possibly sparse), and A is a sparse matrix? Note that this expression is resulting in one single number (scalar). When I try to implement this, it is very very slow...

Do you think there is any way to code this algorithm faster using Eigen structures and available functions? Thanks.