The Jacobi eigenvalue algorithm

This is a classical algorithm proposed by the nineteenth century mathematician C. G. J. Jacobi in connection with some astronomical computations. See wikipedia for a detailed description and some historical references.

The method was computationally tedious, and remained dormant until the advent of modern computers in the mid 20th century. Since its re-discovery it has been refined and improved many times, though much faster algorithms have since been devised and implemented.

I first met the Jacobi algorithm as an early Fortran programming exercise I had as a student in 1966. It’s simplicity and ingenuity fascinated me then and kindled an interest in numerical computations of this kind that has remained ever since. It was a very good way to learn programming.

Parallel revival

There has been some renewed interest in Jacobi-like methods in recent times, however, since unlike the faster methods for eigensolution computations, it offers the possibility of parallelisation. See, for example, Zhou and Brent for one possibility, and others in the references therein.

Purpose of this package

This is a demonstration package used for teaching purposes. It’s main purposes are to provide an example of an intermediate-level programming task where an efficient coding in pure R and one using in C++ using Rcpp are strikingly similar. The task also involves matrix manipulation in pure Rcpp, rather than using RcppArmadillo for example, which is of some teaching interest as well.

There are some situations where the C++ function provided, JacopiCpp, can be slightly faster than the in-built eigen function in the base package, mainly for large numbers of small symmetric matrices. Persons with a fascination for old algorithms might find the comparison with modern versions and alternatives interesting, but generally the functions are not intended for production use.

If someone is motivated to take up the challenge of producing a fast parallel Jacobi algorithm coding in R and provide it as a package, there may well be much practical interest (and this package will have served a useful practical purpose, if somewhat vicariously).

Stagewise protocol

The so-called ‘stagewise’ rotation protocol differs from the standard protocol by traversing the off-diagonal elements of the matrix systematically and rather than simply eliminating the element with maximum absolute value only, it eliminates all elements it encounters in turn if their absolute value exceeds a progressively reducing threshold. This avoids some repeated searching of the off-diagonal section of the matrix and confers a small, if useful, speed advantage. This implementation is provided in JacobiS.

R

This includes the interface to the Rcpp function.

JacobiR <- function(x, symmetric = TRUE, only.values = FALSE,
                    eps = if(!only.values) .Machine$double.eps else
                      sqrt(.Machine$double.eps)) {
  if(!symmetric) 
    stop("only real symmetric matrices are allowed")
  n <- nrow(x)
  H <- if(only.values) NULL else diag(n)
  eps <- max(eps, .Machine$double.eps)

  if(n > 1) {
    lt <- which(lower.tri(x))

    repeat {
      k <- lt[which.max(abs(x[lt]))]  ## the matrix element
      j <- floor(1 + (k - 2)/(n + 1)) ## the column
      i <- k - n * (j - 1)            ## the row

      if(abs(x[i, j]) < eps) break

      Si <- x[, i]
      Sj <- x[, j]

      theta <- 0.5*atan2(2*Si[j], Sj[j] - Si[i])
      c <- cos(theta)
      s <- sin(theta)

      x[i, ] <- x[, i] <- c*Si - s*Sj
      x[j, ] <- x[, j] <- s*Si + c*Sj
      x[i,j] <- x[j,i] <- 0
      x[i,i] <- c^2*Si[i] - 2*s*c*Si[j] + s^2*Sj[j]
      x[j,j] <- s^2*Si[i] + 2*s*c*Si[j] + c^2*Sj[j]
      if(!only.values) {
        Hi <- H[, i]
        H[, i] <- c*Hi - s*H[, j]
        H[, j] <- s*Hi + c*H[, j]
      }
    }
  }
  structure(list(values = as.vector(diag(x)), vectors = H), 
            class = c("Jacobi", "eigen"))
}
##
## The interface function to the Rcpp code
##
Jacobi <- function(x, symmetric = TRUE, only.values = FALSE, eps = 0.0) {
  if(!symmetric) 
    stop("only real symmetric matrices are allowed")
  structure(JacobiCpp(x, only.values, eps), 
            class = c("Jacobi", "eigen"))
}

Rcpp

We begin with one helper function:

#include <Rcpp.h>
using namespace Rcpp;

NumericMatrix Ident(int n) // not exported.
{
    NumericMatrix I(n, n);
    for(int i = 0; i < n; i++) I(i, i) = 1.0;
    return I;
}

// [[Rcpp::export]]
List JacobiCpp(NumericMatrix x, bool only_values = false, double eps = 0.0)
{
    NumericMatrix S(clone(x));
    int nr = S.nrow();
    bool vectors = !only_values;
    NumericMatrix H;

    if(vectors) {
      H = Ident(nr);
    }

    double eps0 =  as<double>((as<List>(Environment::base_env()[".Machine"]))["double.eps"]);
    double tol = eps > eps0 ? eps : eps0;  // i.e. no lower than .Machine$double.eps
    if(only_values & (eps == 0.0)) tol = sqrt(tol); // a lower accuracy is adequate here.

    while(true) {
        double maxS = 0.0;
        int i=0, j=0;
        for(int row = 1; row < nr; row++) {  // find value & position of maximum |off-diagonal|
            for(int col = 0; col < row; col++) {
                double val = fabs(S(row, col));
                if(maxS < val) {
                  maxS = val;
                  i = row;
                  j = col;
                }
           }
        }
        if(maxS <= tol) break;

        NumericVector Si = S(_, i), Sj = S(_, j);

        double theta = 0.5*atan2(2.0*Si(j), Sj(j) - Si(i));
        double s = sin(theta), c = cos(theta);

        S(i, _) = S(_, i) = c*Si - s*Sj;
        S(j, _) = S(_, j) = s*Si + c*Sj;
        S(i, j) = S(j, i) = 0.0;
        S(i, i) = c*c*Si(i) - 2.0*s*c*Si(j) + s*s*Sj(j);
        S(j, j) = s*s*Si(i) + 2.0*s*c*Si(j) + c*c*Sj(j);

      if(vectors) {
            NumericVector Hi = H(_, i);
            H(_, i) = c*Hi - s*H(_, j);
            H(_, j) = s*Hi + c*H(_, j);
      }
    }
    if(vectors) {
      return List::create(_["values"] = diag(S),
                          _["vectors"] = H);
    } else {
      return List::create(_["values"] = diag(S),
                          _["vectors"] = R_NilValue);
    }
}

// Stagewise protocol version
//
// [[Rcpp::export]]
List JacobiSCpp(NumericMatrix x, bool only_values = false, double eps = 0.0)
{
  NumericMatrix S(clone(x));
  int nr = S.nrow();
  bool vectors = !only_values;
  NumericMatrix H;
  
  if(vectors) {
    H = Ident(nr);
  }
  
  double eps0 =  as<double>((as<List>(Environment::base_env()[".Machine"]))["double.eps"]);
  double tol = eps > eps0 ? eps : eps0;  // i.e. no lower than .Machine$double.eps
  if(only_values & (eps == 0.0)) tol = sqrt(tol); // a lower accuracy is adequate here.

  double maxS = 0.0;
  for(int row = 1; row < nr; row++) {  // find value of maximum |off-diagonal|
    for(int col = 0; col < row; col++) {
      double val = fabs(S(row, col));
      maxS = maxS < val ? val : maxS;
    }
  }
  
  while(true) {
    if(maxS <= tol) break;
    double maxS0 = 0.0;
    for(int i = 1; i < nr; i++) {
      for(int j = 0; j < i; j++) {
        double val = fabs(S(i, j));
        maxS0 = maxS0 < val ? val : maxS0;
        if (val > maxS/2.0) {
          NumericVector Si = S(_, i), Sj = S(_, j);
          
          double theta = 0.5*atan2(2.0*Si(j), Sj(j) - Si(i));
          double s = sin(theta), c = cos(theta);
          
          S(i, _) = S(_, i) = c*Si - s*Sj;
          S(j, _) = S(_, j) = s*Si + c*Sj;
          S(i, j) = S(j, i) = 0.0;
          S(i, i) = c*c*Si(i) - 2.0*s*c*Si(j) + s*s*Sj(j);
          S(j, j) = s*s*Si(i) + 2.0*s*c*Si(j) + c*c*Sj(j);
          
          if(vectors) {
            NumericVector Hi = H(_, i);
            H(_, i) = c*Hi - s*H(_, j);
            H(_, j) = s*Hi + c*H(_, j);
          }
        }
      }
    }
    maxS = maxS0;
  }
  if(vectors) {
    return List::create(_["values"] = diag(S),
                        _["vectors"] = H);
  } else {
    return List::create(_["values"] = diag(S),
                        _["vectors"] = R_NilValue);
  }
}