Eigendecomposition times and errors

eigendecomp_parts_bench.cpp

//#include <stan/math.hpp>
#include <stan/math/prim/mat/fun/symmetric_eigensolver.hpp>
#include <chrono>
#include <cstdio>

using namespace std;
using namespace Eigen;
using namespace stan::math;

using Scalar = double;
void eigendecomposition_test(const vector<int>& sizes, int repeats = 1) {
  for (int size : sizes) {
    int multi_run = std::max(1, static_cast<int>(std::pow(500./size,3)));
    for (int i = 0; i < repeats; i++) {
      Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(size, size);
      Matrix<Scalar,Dynamic,Dynamic> packed(size,size);
      Matrix<Scalar,Dynamic,Dynamic> rnd = Matrix<Scalar,Dynamic,Dynamic>::Random(size, size);
      Matrix<Scalar,Dynamic,Dynamic> vecs(size,size);
      Matrix<Scalar,Dynamic,1> vals(size);
      a += a.transpose().eval();
      auto start = std::chrono::steady_clock::now();
      Tridiagonalization<Matrix<Scalar, Dynamic, Dynamic>> tri;
      for(int j=0;j<multi_run;j++) {
        tri.compute(a);
        packed = tri.packedMatrix();
      }
      int eigen_time = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start).count();

      start = std::chrono::steady_clock::now();
      for(int j=0;j<multi_run;j++) {
        stan::math::internal::block_householder_tridiag(a, packed);
      }
      int cpu_time = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start).count();

      Eigen::Matrix<Scalar,Dynamic,1> diagonal = packed.diagonal();
      Eigen::Matrix<Scalar,Dynamic,1> subdiagonal = packed.diagonal(1);
      start = std::chrono::steady_clock::now();
      for(int j=0;j<multi_run;j++) {
        stan::math::internal::tridiagonal_eigensolver(diagonal, subdiagonal, vals, vecs);
      }
      int tri_solve_time = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start).count();

      start = std::chrono::steady_clock::now();
      for(int j=0;j<multi_run;j++) {
        Eigen::internal::computeFromTridiagonal_impl(diagonal, subdiagonal, 30, true, vecs);
      }
      int eigen_tri_solve_time = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start).count();

      start = std::chrono::steady_clock::now();
      for(int j=0;j<multi_run;j++) {
        rnd = tri.matrixQ() /* * rnd*/;
      }
      int eigen_apply_time = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start).count();

      start = std::chrono::steady_clock::now();
      for(int j=0;j<multi_run;j++) {
        stan::math::internal::block_apply_packed_Q(packed, rnd);
      }
      int cpu_apply_time = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start).count();

      double multi_run_dbl = multi_run;
      printf("%d, %lf, %lf, %lf, %lf, %lf, %lf\n", size,
          eigen_time / multi_run_dbl,
          cpu_time / multi_run_dbl,
          eigen_apply_time / multi_run_dbl,
          cpu_apply_time / multi_run_dbl,
          eigen_tri_solve_time / multi_run_dbl,
          tri_solve_time / multi_run_dbl);
    }
  }
}

int main() {
  printf("size, Eigen tri time[ms], My tri time[ms], Eigen apply time[ms], My apply time[ms], Eigen tri solve time[ms], My tri solve time[ms]\n");
  eigendecomposition_test({100,200,500,1000,2000,3000,4000},10);
//  eigendecomposition_test({100,200,500,1000},3);
//  eigendecomposition_test({4000},5);
}

eigendecomposition_speed_precision.cpp

//#include <stan/math.hpp>
#include <stan/math/prim/mat/fun/symmetric_eigensolver.hpp>
#include <chrono>
#include <cstdio>

using namespace std;
using namespace Eigen;
using namespace stan::math;

using Scalar = double;

double residual(const Matrix<Scalar,Dynamic,Dynamic>& a, const Matrix<Scalar,Dynamic,1>& eigenvals, const Matrix<Scalar,Dynamic,Dynamic>& eigenvecs) {
  return (a * eigenvecs - eigenvecs * eigenvals.asDiagonal()).array().abs().colwise().sum().maxCoeff();
}

double ortho_loss(const Matrix<Scalar,Dynamic,Dynamic>& eigenvecs) {
  return (eigenvecs.transpose() * eigenvecs - Matrix<Scalar,Dynamic,Dynamic>::Identity(eigenvecs.rows(), eigenvecs.rows())).array().abs().maxCoeff();
}

void eigendecomposition_test(const vector<int>& sizes, int repeats = 1) {
  for (int size : sizes) {
    for (int i = 0; i < repeats; i++) {
      Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(size, size);
      Matrix<Scalar,Dynamic,Dynamic> vecs(size,size);
      Matrix<Scalar,Dynamic,1> vals(size);
      a += a.transpose().eval();
      auto start = std::chrono::steady_clock::now();
      SelfAdjointEigenSolver <Matrix<Scalar,Dynamic,Dynamic>> slv(a);
      vals = slv.eigenvalues();
      vecs = slv.eigenvectors();
      int eigen_time = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start).count();
      double eigen_residual = residual(a, vals, vecs);
      double eigen_ortho_loss = ortho_loss(vecs);

      start = std::chrono::steady_clock::now();
      selfadjoint_eigensolver(a, vals, vecs);
      int cpu_time = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::steady_clock::now() - start).count();
      double cpu_residual = residual(a, vals, vecs);
      double cpu_ortho_loss = ortho_loss(vecs);

      printf("%d, %d, %d, %e, %e, %e, %e\n", size, eigen_time, cpu_time, eigen_residual, cpu_residual, eigen_ortho_loss, cpu_ortho_loss);
    }
  }
}

int main() {
  printf("tri, tri_eig, size, Eigen time[ms], My time[ms], Eigen max residual, My max residual, Eigen orthogonality loss, My orthogonality loss\n");
  eigendecomposition_test({100,200,500,1000,2000,3000,4000},10);
//  eigendecomposition_test({4000},5);
}

eigendecomposition_speed_precision.csv

parts_times.csv

plotRes.py

import pandas, numpy
import matplotlib.pyplot as plt


data = pandas.read_csv("parts_times.csv",sep = ", ")
size = numpy.array(data["size"])

data2 = pandas.read_csv("eigendecomposition_speed_precision.csv",sep = ", ")
size2 = data2["size"]

plt.plot([0, size[-1]], [1,1])
plt.plot(size, data["Eigen tri time[ms]"] / data["My tri time[ms]"], "o", label="tridiagonalization")
plt.plot(size, data["Eigen apply time[ms]"] / data["My apply time[ms]"], "o", label="Multiply w Q / coonstruct Q")
plt.plot(size, data["Eigen tri solve time[ms]"] / data["My tri solve time[ms]"], "o", label="MRRR / QR iter")

plt.plot(size2, data2["Eigen time[ms]"] / data2["My time[ms]"], "o", label="whole algorithm combined")
plt.xlabel("size")
plt.ylabel("speedup")
plt.title("Speedup")
plt.legend()
plt.show()


for col in ("Eigen max residual", "My max residual", "Eigen orthogonality loss", "My orthogonality loss"):
    plt.plot(size2, data2[col], "o", label=col)
plt.xlabel("size")
plt.ylabel("error")
plt.title("Errors")
plt.semilogy()
plt.legend()
plt.show()