eigenvalue_estimation.cpp

#include <cassert>
#include <cmath>
#include <fstream>
#include <sstream>
#include <iostream>
#include <random>
#include <string>
#include <vector>
#include <chrono>
#include <iomanip>

// ---------------- Utilities ----------------
void output_error_msg(const std::string& msg)
{
    std::cerr << "Error: " << msg << std::endl;
    std::exit(EXIT_FAILURE);
}

// ---------------- Linear Algebra Ops ----------------
// Dot Product: x^T * y
double dot(const std::vector<double>& a,
           const std::vector<double>& b)
{
    assert(a.size() == b.size());
    double s = 0.0;
    for (size_t i = 0; i < a.size(); ++i)
        s += a[i] * b[i];
    return s;
}

double norm2(const std::vector<double>& x)
{
    return std::sqrt(dot(x, x));
}

void normalize(std::vector<double>& x)
{
    double n = norm2(x);
    if (n == 0.0)
        output_error_msg("Zero vector encountered during normalization.");
    for (double& v : x)
        v /= n;
}

// Vector Update (AXPY): y = a * x + y
void axpy(double alpha,
          const std::vector<double>& x,
          std::vector<double>& y)
{
    assert(x.size() == y.size());
    for (size_t i = 0; i < x.size(); ++i)
        y[i] += alpha * x[i];
}

std::vector<double> matvec(const std::vector<double>& A,
                           const std::vector<double>& x,
                           int n)
{
    assert((int)x.size() == n);
    std::vector<double> y(n, 0.0);
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            y[i] += A[i * n + j] * x[j];
    return y;
}

// ---------------- Log data during iterating ----------------
struct IterationLogger
{
    std::ofstream file;

    explicit IterationLogger(const std::string& filename)
    {
        file.open(filename);
        if (!file)
            output_error_msg("Failed to open log file " + filename);
        file << "Iter, Eigenvalue, Relative_residual \n";
    }

    void log(size_t iter, double eig, double res)
    {
        file << iter << "," << eig << "," << res << "\n";
    }
};

// ---------------- Preconditioned Conjugate Gradient Solver ----------------
int pcg(const std::vector<double>& A,
        const std::vector<double>& b,
        std::vector<double>& x,
        int n,
        const int cg_max_iter,
        const double cg_tol)
{
    std::vector<double> r(n, 0.0), z(n, 0.0), p(n, 0.0), Ap(n, 0.0);
    std::vector<double> Minv(n, 0.);

    // Jacobi preconditioner: M^-1 = 1/diag(A)
    for (int i = 0; i < n; ++i)
        Minv[i] = 1.0 / A[i * n + i];

    Ap = matvec(A, x, n);

    for (int i = 0; i < n; ++i)
        r[i] = b[i] - Ap[i];

    const double bnorm = norm2(b);

    // Preconditioning: z = M^-1 * r
    for (int i = 0; i < n; ++i)
        z[i] = Minv[i] * r[i];

    p = z;
    double rz_old = dot(r, z);

    for (int k = 0; k < cg_max_iter; ++k)
    {
        Ap = matvec(A, p, n);

        // check for divison by ~ 0
        if (std::abs(dot(p, Ap)) < 1e-30)
            break;

        const double alpha = rz_old / dot(p, Ap);

        axpy(alpha, p, x); // x = x + alpha*p
        axpy(-alpha, Ap, r); // r = r - alpha*Ap

        const double rnorm = norm2(r);

        if (rnorm / bnorm < cg_tol)
            return k + 1;

        // Apply Preconditioner
        for (int i = 0; i < n; ++i)
            z[i] = Minv[i] * r[i];

        const double rz_new = dot(r, z);

        // check for divison by ~ 0
        if (std::abs(rz_old) < 1e-30)
            break;

        const double beta = rz_new / rz_old;

        // p = z + beta*p
        for (int i = 0; i < n; ++i)
            p[i] = z[i] + beta * p[i];

        rz_old = rz_new;
    }
    return cg_max_iter;
}

// ---------------- Eigenvalue Estimators ----------------
// lambda_max ≈ rho = (x^T * A * x) / (x^T * x)
double rayleigh_quotient(const std::vector<double>& x,
                         const std::vector<double>& Ax,
                         int n)
{
    return dot(x, Ax) / dot(x, x);
}

// ---------------- Relative residuals in computing eigenvalue ----------------
double relative_eigen_residual(const std::vector<double>& A,
                               const std::vector<double>& x,
                               int n,
                               double lambda)
{
    std::vector<double> Ax = matvec(A, x, n);

    double num = 0.0, denom = 0.0;

    // Residual r = A*x - lambda*x
    for (int i = 0; i < n; ++i)
    {
        double r = Ax[i] - lambda * x[i];

        // Relative residual = norm(r) / norm(Ax)
        num += r * r;
        denom += Ax[i] * Ax[i];
    }

    return std::sqrt(num / denom);
}

// ---------------- Power iteration to estimate lambda_max ----------------
double power_method(const std::vector<double>& A,
                    std::vector<double>& x,
                    int n,
                    const int maxit,
                    const double tol,
                    IterationLogger& log)
{
    // Start timer
    auto t_start = std::chrono::steady_clock::now();
    double solve_time = 0.0;

    normalize(x);

    std::vector<double> Ax(x.size(), 0.0);
    double lambda_old = 0.0;

    for (int k = 0; k < maxit; ++k)
    {
        Ax = matvec(A, x, n);
        double lambda = rayleigh_quotient(x, Ax, n);

        // Compute relative residual in eigenvalue
        double res = relative_eigen_residual(A, x, n, lambda);

        log.log(k, lambda, res);

        normalize(Ax);
        x = Ax;

        if (res < tol)
        {
            std::cout << "\nPower method converged in " << k+1 << " iterations\n";
            auto t_end = std::chrono::steady_clock::now();
            solve_time = std::chrono::duration<double>(t_end - t_start).count();
            std::cout << "Solve time for power method: " << solve_time << " seconds\n";
            return lambda;
        }

        if (k == maxit - 1)
            std::cerr << "Warning: Max iterations reached in power method. Residual = " << res << "\n";

        lambda_old = lambda;
    }

    auto t_end = std::chrono::steady_clock::now();
    solve_time = std::chrono::duration<double>(t_end - t_start).count();
    std::cout << "Solve time for power method: " << solve_time << " seconds\n";
    return lambda_old;
}

// ------- Inverse power iteration to estimate lambda_min via PCG solver -------
double inverse_iteration(const std::vector<double>& A,
                         std::vector<double>& x,
                         int n,
                         const int cg_max_iter,
                         const double cg_tol,
                         const int maxit,
                         const double tol,
                         IterationLogger& log)
{
    // Start timer
    auto t_start = std::chrono::steady_clock::now();
    double solve_time = 0.0;

    normalize(x);

    std::vector<double> Ax(x.size(), 0.0);
    double lambda_old = 0.0;

    for (int k = 0; k < maxit; ++k)
    {
        // Solve A y = x (approximately)
        std::vector<double> y(n, 0.0);

        // Use PCG to approximate y_k ~ Ainv * x_k
        // with Jacobi preconditioner Minv = 1/diag(A)
        int converged_cg_iter_num = pcg(A, x, y, n, cg_max_iter, cg_tol);

        if (converged_cg_iter_num == cg_max_iter)
            std::cerr << "Warning: CG did not converge at iteration " << k << "\n";

        normalize(y);
        x = y;

        Ax = matvec(A, x, n);

        double lambda = rayleigh_quotient(x, Ax, n);

        // Compute relative residual in eigenvalue
        double res = relative_eigen_residual(A, x, n, lambda);

        log.log(k, lambda, res);

        if (res < tol)
        {
            std::cout << "\nInverse power method converged in " << k+1 << " iterations\n";
            auto t_end = std::chrono::steady_clock::now();
            solve_time = std::chrono::duration<double>(t_end - t_start).count();
            std::cout << "Solve time for inverse power method: " << solve_time << " seconds\n";
            return lambda;
        }

        if (k == maxit - 1)
            std::cerr << "Warning: Max iterations reached in inverse power method. Residual = " << res << "\n";

        lambda_old = lambda;
    }

    auto t_end = std::chrono::steady_clock::now();
    solve_time = std::chrono::duration<double>(t_end - t_start).count();
    std::cout << "Solve time for inverse power method: " << solve_time << " seconds\n";
    return lambda_old;
}

// ---------------- Gershgorin Bounds ----------------
// lamba \in [ Aii - Ri, Aii + Ri ]
// lambda_min >= min_i( Aii - Ri ), lambda_max <= max_i(Aii + Ri )
void gershgorin_bounds(const std::vector<double>& A,
                       int n,
                       double& gmin,
                       double& gmax)
{
    gmin = std::numeric_limits<double>::max();
    gmax = 0.0;

    for (int i = 0; i < n; ++i)
    {
        double center = A[i * n + i]; // Diagonal entry Aii
        double radius = 0.0;

        for (int j = 0; j < n; ++j)
            if (j != i)
                radius += std::abs(A[i * n + j]); // Ri = sum_(j!=i) abs(Aij)

        gmin = std::min(gmin, center - radius);
        gmax = std::max(gmax, center + radius);
    }
}

// --------- Helper to load input CSV files -------------
std::vector<double> read_matrix_csv(const std::string& file,
                                    int& n)
{
    std::ifstream in(file);
    if (!in)
        output_error_msg("Cannot open matrix file: " + file);

    std::vector<double> Matrix;
    std::string line;

    while (std::getline(in, line))
    {
        size_t pos = 0;
        while ((pos = line.find(',')) != std::string::npos)
        {
            Matrix.push_back(std::stod(line.substr(0, pos)));
            line.erase(0, pos + 1);
        }
        Matrix.push_back(std::stod(line));
    }

    n = std::sqrt(Matrix.size());
    if (n * n != (int)Matrix.size())
        output_error_msg("Matrix is not square.");

    return Matrix;
}

std::vector<double> read_vector_csv(const std::string& filename)
{
    std::ifstream in(filename);
    if (!in)
        output_error_msg("Cannot open vector file: " + filename);

    std::vector<double> vec;
    std::string line;

    while (std::getline(in, line))
    {
        if (line.empty()) continue;

        std::stringstream ss(line);
        double value;
        ss >> value;

        if (ss.fail())
            output_error_msg("Invalid numeric entry in vector file: " + filename);

        vec.push_back(value);
    }

    if (vec.empty())
        output_error_msg("Vector file is empty: " + filename);

    return vec;
}

// ---------------- Eigenvalue Verification using LAPACK ---------------
// LAPACK function declaration (Fortran interface)
// Computes all eigenvalues and optionally eigenvectors of a real symmetric matrix
extern "C" {
void dsyev_(
    char* jobz,      // 'N' = eigenvalues only, 'V' = eigenvalues + eigenvectors
    char* uplo,      // 'U' = upper triangle, 'L' = lower triangle
    int* n,          // Order of matrix
    double* A,       // Matrix (overwritten on output)
    int* lda,        // Leading dimension of A
    double* W,       // Eigenvalues in ascending order
    double* work,    // Workspace
    int* lwork,      // Size of work array
    int* info        // Status: 0=success, <0=illegal arg, >0=convergence failure
    );
}

// Use only for small matrix size O(n^3)
void verify_eigenvalues_lapack(const std::vector<double>& A,
                               int n,
                               double lambda_max_computed,
                               double lambda_min_computed,
                               double kappa_computed)
{
    std::cout << "\nVerification: Computing exact eigenvalues using LAPACK\n";

    // Make a copy of A (LAPACK destroys input matrix)
    std::vector<double> A_copy = A;

    // Eigenvalue storage
    std::vector<double> eigenvalues(n);

    // LAPACK parameters
    char jobz = 'N';  // Only compute eigenvalues, not eigenvectors
    char uplo = 'U';  // Use upper triangle
    int lda = n;
    int info;

    // Query optimal workspace size
    int lwork = -1;
    double work_query;
    dsyev_(&jobz, &uplo, &n, A_copy.data(), &lda, eigenvalues.data(),
           &work_query, &lwork, &info);

    // Allocate workspace
    lwork = (int)work_query;
    std::vector<double> work(lwork);

    // Compute eigenvalues
    dsyev_(&jobz, &uplo, &n, A_copy.data(), &lda, eigenvalues.data(),
           work.data(), &lwork, &info);

    if (info != 0)
    {
        std::cerr << "LAPACK dsyev failed with info = " << info << "\n";
        return;
    }

    // Eigenvalues are returned in ascending order
    double lambda_min_exact = eigenvalues[0];
    double lambda_max_exact = eigenvalues[n - 1];
    double kappa_exact = lambda_max_exact / lambda_min_exact;

    // Compute relative errors
    double error_max = std::abs(lambda_max_computed - lambda_max_exact) / lambda_max_exact;
    double error_min = std::abs(lambda_min_computed - lambda_min_exact) / lambda_min_exact;
    double error_kappa = std::abs(kappa_computed - kappa_exact) / kappa_exact;

    // Display results
    std::cout << "\nExact values (LAPACK):\n";
    std::cout << "  lambda_max (exact) = " << std::setprecision(10) << lambda_max_exact << "\n";
    std::cout << "  lambda_min (exact) = " << std::setprecision(10) << lambda_min_exact << "\n";
    std::cout << "  kappa (exact) = " << std::setprecision(10) << kappa_exact << "\n";

    std::cout << "\nRelative errors:\n";
    std::cout << "  Relative error in lambda_max = " << std::fixed << error_max * 100 << "%\n";
    std::cout << "  Relative error in lambda_min = " << std::fixed << error_min * 100 << "%\n";
    std::cout << "  Relative error in kappa = " << std::fixed << error_kappa * 100 << "%\n";

    // Validation checks
    std::cout << "\nValidation:\n";
    if (error_max < 1e-6 && error_min < 1e-6)
    {
        std::cout << "Eigenvalues match to 6+ digits.\n";
    }
    else if (error_max < 1e-4 && error_min < 1e-4)
    {
        std::cout << "Eigenvalues match to 4+ digits.\n";
    }
    else
    {
        std::cout << "Eigenvalue errors exceed expected tolerance.\n";
    }
}

// ---------------- Main ----------------
int main(int argc, char* argv[])
{
    if (argc != 3)
        output_error_msg("Usage: ./eigenvalue_estimation A.csv b.csv");

    int n;

    // Read Matrix A and RHS vector b
    std::vector<double> A = read_matrix_csv(argv[1], n);
    std::vector<double> b = read_vector_csv(argv[2]);

    assert(n > 0);
    assert(b.size() == n && "Dimension mismatch: A and b");
    std::cout << "Matrix size: " << n << "x" << n << "\n";

    // Compute Gershgorin bounds
    double gmin, gmax;
    gershgorin_bounds(A, n, gmin, gmax);

    std::cout << "\nGershgorin bounds: [" << gmin << ", " << gmax << "]\n";

    // ---- Initial random vectors ----
    //std::random_device rd;  // Used to obtain a seed for random number engine
    // Standard mersenne_twister_engine seeded with fixed seed
    // alternately can use rd() as random seed
    std::mt19937 gen(42);
    std::uniform_real_distribution<double> dist(0.0, 1.0);

    std::vector<double> x_max(n), x_min(n);
    for (int i = 0; i < n; ++i)
    {
        x_max[i] = dist(gen);
        x_min[i] = dist(gen);
    }

    if (norm2(x_max) < 1e-14 || norm2(x_min) < 1e-14)
        output_error_msg("Zero initial vector generated.");

    // Iteration logger for power and inverse power methods
    IterationLogger log_max("lambda_max.csv");
    IterationLogger log_min("lambda_min.csv");

    int max_iter = 500;
    double tol = 1e-8;
    int cg_max_iter = 1000;
    double cg_tol = 0.1 * tol; // Need high CG tol for large condition number cases

    double lambda_max = power_method(A, x_max, n, max_iter, tol, log_max);
    double lambda_min = inverse_iteration(A, x_min, n, cg_max_iter, cg_tol, max_iter, tol, log_min);
    double kappa = lambda_max / lambda_min;

    // Sanity checks for positive eigen values
    if (lambda_min <= 0 || lambda_max <= 0)
        output_error_msg("Non-positive eigenvalue detected. Matrix may not be SPD.");

    // Sanity check for eigen values within Gershgorin bounds
    if(lambda_min < gmin)
        std::cerr <<"Minimum eigenvalue lambda_min below the minimum Gershgorin bound.\n";
    if(lambda_max > gmax)
        std::cerr << "Maximum eigenvalue lambda_max above the maximum Gershgorin bound.\n";

    std::cout << "\nEstimated values:\n";
    std::cout << "  Maximum eigenvalue lambda_max ≈ " << std::setprecision(10) << lambda_max << "\n";
    std::cout << "  Minimum eigenvalue lambda_min ≈ " << std::setprecision(10) << lambda_min << "\n";
    std::cout << "  Condition number kappa ≈ " << std::setprecision(10) << kappa << "\n";

    // Verify against LAPACK (only for small matrices)
    if (n <= 64)
    {
        verify_eigenvalues_lapack(A, n, lambda_max, lambda_min, kappa);
    }
    else
    {
        std::cout << "\nSkipping LAPACK verification: matrix too large.\n";
    }

    return 0;
}