Simple approximate eigensolver for paniq

gistfile1.cpp

static float const singularEps = 2 * 1.192092896e-7f; // float32 epsilon, *2 to have a bit of margin

// QR factorization using MGS
// store 1 / diag elements in diag of R since that's what we need
bool QR(IN(mat3,m),OUT(mat3,q),OUT(mat3,r)) {
    q[0] = normalize(m[0]);
    q[1] = normalize(m[1] - dot(m[1], q[0])*q[0]);
    q[2] = cross(q[0], q[1]);

    float d0 = dot(m[0], q[0]);
    float d1 = dot(m[1], q[1]);
    float d2 = dot(m[2], q[2]);
    float maxd = max(max(fabs(d0), fabs(d1)), fabs(d2));
    float mind = min(min(fabs(d0), fabs(d1)), fabs(d2));

    // Are we numerically singular? (This test is written to work
    // in the presence of NaN/Inf; using >= won't work right.)
    if (!(maxd * singularEps < mind))
        return false;
    
    r[0] = vec3(1.0 / d0, 0.0, 0.0);
    r[1] = vec3(dot(m[1],q[0]), 1.0 / d1, 0.0);
    r[2] = vec3(dot(m[2],q[0]), dot(m[2],q[1]), 1.0 / d2);
    return true;
}

// matrix a must be upper triangular
// with main diagonal storing reciprocal of actual vals
vec3 solve_Ux_b(IN(mat3,a), IN(vec3,b)) {
    float x2 = b[2] * a[2][2];
    float x1 = (b[1] - a[2][1]*x2) * a[1][1];
    float x0 = (b[0] - a[1][0]*x1 - a[2][0]*x2) * a[0][0];
    return vec3(x0,x1,x2);
}    

// rayleigh quotient iteration from Q R matrices
void rayleigh_quot(IN(mat3,a), INOUT(float,mu), INOUT(vec3, x)) {
    mat3 q, r;
    vec3 y = x;
    for (int i = 0; i < 5; ++i) {
        x = y / length(y);
        if (!QR(a - mat3(mu), q, r))
            break;
        y = solve_Ux_b(r, x * q);
        mu = mu + 1.0 / dot(y,x);
    }
}

// inverse iter to find EV closest to mu
void inverse_iter(IN(mat3,a), INOUT(float,mu), INOUT(vec3, x)) {
    mat3 q, r;

    // If a - mat3(mu) is singular, we already have an eigenvector!
    if (!QR(a - mat3(mu), q, r))
        return;

    // If you know eigenvalues aren't too large/small, can skip
    // normalize here. (It's only there to present over-/underflow.)
    // Normalizing once at the end (before calc of mu) works fine.
    for (int i = 0; i < 4; ++i)
        x = normalize(solve_Ux_b(r, x * q));

    mu = dot(x, a*x);
}

// ---- test

void main() {
    mat3 a = mat3(
        1.0,1.0,3.0,
        2.0,2.0,2.0,
        3.0,1.0,1.0);

    vec3 x = vec3(1.0);
    float mu = 0;

    // so the problem with rayleigh is that it really likes to latch on to
    // the eigenvalue with largest absolute value.
    // better to use normal inverse iteration if we want the EV closest to 0!
    //rayleigh_quot(a, mu, x);

    // once to get close to smallest EV, second pass to polish
    inverse_iter(a, mu, x);
    inverse_iter(a, mu, x);
    
    printf("mu=%.5f\n",mu);
    printf("x="); dump_vec3(x);
}