Skip to content

Iterative — Implementation

File Structure

iterative/
├── tiny_iterative.h        (35 lines — API declarations)
└── tiny_iterative.c        (755 lines — implementation)

Dependencies: tiny_eigen.h, tiny_decomp.h, tiny_constants.h, <math.h>.

Implementation Summary

Component Lines Key Algorithm
Operator init helpers 80 Dense / Hankel / Block-Hankel / Covariance
Arnoldi 120 MGS-based Krylov iteration
Lanczos 80 Tridiagonal symmetric iteration
RSVD core 250 Random projection + QR + small SVD
SVD utilities 20 Energy / rank / condition
Random number generation 50 Box-Muller Gaussian RNG

Operator Interface

Each operator init function creates a tiny_operator_f32_t with matvec and rmatvec callbacks:

// Dense operator: matvec = A*x, rmatvec = A^T*x
tiny_operator_f32_t tiny_dense_operator_init(const float *A, int m, int n, int step)
{
    tiny_operator_f32_t op = {m, n, NULL, NULL, NULL};
    op.ctx = (void*)&((struct dense_ctx){A, m, n, step});
    op.matvec = dense_matvec;
    op.rmatvec = dense_rmatvec;
    return op;
}

Hankel operator builds an implicit Hankel matrix from a 1D signal. Instead of constructing the full \(m \times n\) matrix (\(O(mn)\) memory), matvec computes \(y_k = \sum_{j=0}^{n-1} s_{k+j} \cdot x_j\) on-the-fly (\(O(m \cdot n)\) time, \(O(1)\) extra memory).

Arnoldi: Core Loop

for (int j = 0; j < k; j++) {
    // 1. Apply operator to current basis vector
    op->matvec(&V[(j)*n], &V[(j+1)*n], op->ctx);   // v_{j+1} = A * v_j

    // 2. MGS orthogonalize against all previous vectors
    for (int i = 0; i <= j; i++) {
        float h = dot_product(&V[(j+1)*n], &V[i*n], n);  // h_{i,j} = v_i · v_{j+1}
        H[i*k + j] = h;
        axpy(&V[(j+1)*n], -h, &V[i*n], n);              // v_{j+1} -= h * v_i
    }

    // 3. Normalize
    float norm = sqrt(dot_product(&V[(j+1)*n], &V[(j+1)*n], n));
    H[(j+1)*k + j] = norm;
    if (norm > 0) scale(&V[(j+1)*n], 1.0f / norm, n);
    else break;  // invariant subspace found
}

The Hessenberg matrix \(H\) is reused from subsequent eigenvalue computation (e.g., finding Ritz values).

Lanczos: Tridiagonal Variant

Lanczos exploits symmetry (\(A = A^T\)), producing only \(\alpha\) (diagonal) and \(\beta\) (off-diagonal):

for (int j = 0; j < k; j++) {
    // 1. Apply operator
    op->matvec(&V[j*n], w, op->ctx);

    // 2. α_{j} = v_j · (A·v_j)
    alpha[j] = dot_product(&V[j*n], w, n);

    // 3. Orthogonalize against v_{j-1} and v_j
    if (j > 0) axpy(w, -beta[j-1], &V[(j-1)*n], n);
    axpy(w, -alpha[j], &V[j*n], n);

    // 4. β_j = ||w||
    beta[j] = sqrt(dot_product(w, w, n));
    if (beta[j] > 0) copy(&V[(j+1)*n], w, n); scale(&V[(j+1)*n], 1/beta[j], n);
    else break;
}

Only 2 storage vectors needed vs Arnoldi's \(k\) vectors — significant savings for large \(n\).

Randomized SVD

Input:  A (m×n), target rank k, oversampling p, power q
Output: U (m×k), S (k), Vt (k×n)
──────────────────────────────────────────────
1. Ω = randn(n, k+p)              // Gaussian random matrix
2. Y = A · Ω                      // m × (k+p)
3. Y = QR → Q                     // thin QR
4. B = Q^T · A                    // (k+p) × n
5. [Ũ, Σ, V] = svd(B)            // full SVD of tiny B
6. U = Q · Ũ                      // reconstruct left singular vectors

Power iteration (step 4 in RSVD):

for (int i = 0; i < q; i++) {
    // Y = A * (A^T * Y) — power iteration amplifies dominant singular values
    op->rmatvec(Y, temp, op->ctx);   // temp = A^T * Y
    op->matvec(temp, Y, op->ctx);    // Y = A * temp
}

Each \(q\)-pass improves the approximation for singular values that decay slowly.

SVD Utilities

All utilities operate on the singular value array \(S\):

// Cumulative energy: energy[i] = sum_{j=0}^{i} S_j^2 / sum(S^2)
void tiny_svd_cumulative_energy_f32(const float *S, int len, float *energy) {
    float total = 0;
    for (int i = 0; i < len; i++) total += S[i]*S[i];
    float acc = 0;
    for (int i = 0; i < len; i++) {
        acc += S[i]*S[i];
        energy[i] = acc / total;
    }
}

// Model order selection
int order = tiny_svd_rank_by_energy_f32(S, len, 0.99f);  // 99% energy
int rank  = tiny_svd_rank_by_threshold_f32(S, len, 1e-3f); // numerical rank
float kappa = tiny_svd_condition_estimate_f32(S, len);      // condition number