Matrix Operations / Sparse Matrices

Why Matrix Operations / Sparse Matrices Matters

A network model may contain millions of variables while each equation involves only a few neighbors. Its matrix has mostly zero entries. Storing every zero wastes memory, and multiplying through them wastes work. Sparse representations store the meaningful entries and their locations, making operations such as matrix–vector multiplication practical.

Solving a system introduces another concern: numerical stability. A zero diagonal entry does not necessarily mean the system is singular; exchanging rows can provide a usable pivot. Elimination can also create new nonzeros, called fill-in, so a sparse solver cannot assume that the initial nonzero pattern stays fixed.

This lesson connects basic matrix operations, sparse storage, sparse matrix–vector products, and a reference linear-system solver. You will trace Gaussian elimination with partial pivoting, implement it using plain nested arrays, and check its answer against the original equations. A separate row-list routine demonstrates sparse multiplication. The solver deliberately uses dense workspace so every arithmetic step is visible; the sparse-storage discussion explains what changes at larger scale.

Core Concepts and Terminology

An $m\times n$ matrix has $m$ rows and $n$ columns. Addition requires equal shapes, with $(A+B)_{ij}=A_{ij}+B_{ij}$. Transposition exchanges the indices: $(A^{\mathsf T})_{ji}=A_{ij}$. If $A$ is $m\times n$ and $B$ is $n\times q$, their product is $m\times q$ with $(AB)_{ij}=\sum_{k=0}^{n-1}A_{ik}B_{kj}$.

For a vector $x$ of length $n$, matrix–vector multiplication computes $y_i=\sum_j A_{ij}x_j$. Solving $Ax=b$ asks for an unknown $x$; it is a different operation from multiplying a known vector. Our solve routine accepts a nonempty square matrix and a matching right-hand side.

A sparse matrix has sufficiently few stored entries for sparse storage to be useful. Let $s=\operatorname{nnz}(A)$ be the number of nonzero entries. Common formats include coordinate triples, or COO; compressed sparse row, or CSR; and row lists. CSR uses a value array, column indices, and row pointers, as described in Netlib's compressed-row storage reference. Row lists store each row's column indices and values separately; the programs use this simple nested-array form for multiplication.

Use zero-based indices throughout. The worked matrix is

$$ A=\begin{pmatrix}0&2&1\\1&3&0\\4&0&5\end{pmatrix},\qquad b=\begin{pmatrix}7\\7\\19\end{pmatrix}. $$

Its six nonzeros have row-list columns [[1, 2], [0, 1], [0, 2]] and corresponding values [[2, 1], [1, 3], [4, 5]]. The equivalent zero-based CSR arrays are row_ptr = [0, 2, 4, 6], col_idx = [1, 2, 0, 1, 0, 2], and values = [2, 1, 1, 3, 4, 5]. Empty rows have empty row lists or equal consecutive CSR pointers.

How Matrix Operations / Sparse Matrices Works, Step by Step

The sparse multiplication routine visits only the listed entries in each row. Start its output at zero, then add the stored value times the vector entry at the stored column. Within a row, the accumulator always equals the contribution of entries already visited. This requires matching column/value row lengths and valid column indices.

The dense reference solver uses partial pivoting, the same pivoting principle used by general LU-based solvers such as LAPACK's DGESV. Its steps are:

1. Validate a nonempty square shape, matching right-hand-side length, finite input values, and a finite tolerance $0\lt\tau\lt1$. Copy $A$ into workspace $U$ and $b$ into a working vector $g$.
2. Set the scale to $\alpha=\max_{ij}\lvert A_{ij}\rvert$ and the threshold to $\tau\alpha$. The default relative tolerance is $\tau=10^{-12}$.
3. For column $k$, find the first row attaining the largest $\lvert U_{ik}\rvert$ among $i\ge k$. Reject if that pivot magnitude is at most the threshold.
4. Swap the selected row into position $k$ in both $U$ and $g$.
5. For every lower row $i$, compute $f=U_{ik}/U_{kk}$. Set $U_{ik}=0$ explicitly, update each later column by $U_{ij}\leftarrow U_{ij}-fU_{kj}$, and update $g_i\leftarrow g_i-fg_k$.
6. After all columns, solve the upper-triangular system from the last row upward using

$$ x_i=\frac{g_i-\sum_{j=i+1}^{n-1}U_{ij}x_j}{U_{ii}}. $$

1. Evaluate the residual with the original $A$ and $b$. A residual computed from overwritten elimination workspace would not independently check the original equations.

The threshold is a numerical acceptance rule, not an exact rank test. It can reject an invertible but poorly scaled system. Partial pivoting controls the current multipliers, but neither pivoting nor this tolerance guarantees an accurate answer for every ill-conditioned problem.

flowchart TD
    A["Read A, b, and tolerance"] --> B{"Valid finite square system?"}
    B -- No --> C["Reject input"]
    B -- Yes --> D["Copy system; compute scale; set k = 0"]
    D --> E{"k less than n?"}
    E -- Yes --> F["Choose largest absolute active-column pivot"]
    F --> G{"Pivot exceeds threshold?"}
    G -- No --> H["Report singular or numerically rejected system"]
    G -- Yes --> I["Swap matrix rows and right-hand-side entries"]
    I --> J["Eliminate every lower row, retaining fill-in"]
    J --> K["Increment k"]
    K --> E
    E -- No --> L["Back substitute"]
    L --> M["Check original-system residual and sparse product"]
    M --> N["Return solution and diagnostics"]

Worked Example with a Full Trace

Start with

$$ [A\mid b]=\left[\begin{array}{ccc|c}0&2&1&7\\1&3&0&7\\4&0&5&19\end{array}\right]. $$

The scale is $\alpha=5$ and the threshold is $5\times10^{-12}$. At $k=0$, the pivot candidates have magnitudes $0,1,4$. Select row $2$ and swap it with row $0$, including the right-hand side:

$$ [U\mid g]=\left[\begin{array}{ccc|c}4&0&5&19\\1&3&0&7\\0&2&1&7\end{array}\right]. $$

For row $1$, the factor is $1/4$. The later-column updates are $3-(1/4)0=3$ and $0-(1/4)5=-5/4$; the right-hand side becomes $7-(1/4)19=9/4$. This creates fill-in in column $2$. Row $2$ has factor $0$ and stays unchanged. After this column,

$$ [U\mid g]=\left[\begin{array}{ccc|c}4&0&5&19\\0&3&-5/4&9/4\\0&2&1&7\end{array}\right]. $$

At $k=1$, compare magnitudes $3$ and $2$. Keep row $1$, with pivot $3$. The only lower factor is $2/3$. Set its column-one entry to zero, update its final entry to $1-(2/3)(-5/4)=11/6$, and update its right-hand side to $7-(2/3)(9/4)=11/2$:

$$ [U\mid g]=\left[\begin{array}{ccc|c}4&0&5&19\\0&3&-5/4&9/4\\0&0&11/6&11/2\end{array}\right]. $$

At $k=2$, the last pivot is $11/6$, above the threshold, and there are no lower rows. Back substitution gives $x_2=(11/2)/(11/6)=3$, then $x_1=(9/4-(-5/4)3)/3=2$, then $x_0=(19-0\cdot2-5\cdot3)/4=1$.

Check the sparse product using the original row lists. Row $0$ accumulates $2\times2=4$, then $4+1\times3=7$. Row $1$ accumulates $1\times1=1$, then $1+3\times2=7$. Row $2$ accumulates $4\times1=4$, then $4+5\times3=19$. Thus $Ax=[7,7,19]^{\mathsf T}=b$, and this sample's computed residual norm is zero. General floating-point problems will usually have small nonzero residuals.

Implementation in Python, C++, and Java

Each panel contains a pivoted dense solver, sparse row-list multiplication, and a residual check. The solver returns a fresh solution and preserves its inputs. The sparse helper treats each row as corresponding column and value lists; duplicate column entries add together, and empty rows are allowed. Its inputs must contain finite vector and stored values. This helper is a sparse multiplication kernel; the elimination solver uses dense nested arrays and retains all fill-in.

All programs use double-precision arithmetic and tolerance $10^{-12}$. They print the same five lines: nnz = 6; x = [1.000000, 2.000000, 3.000000]; sparse Ax = [7.000000, 7.000000, 19.000000]; residual_inf = 0.000e+00; residual check = PASS. The residual check uses an absolute threshold of $10^{-10}$ for this small sample. Save the Java panel as Main.java.

For a general application, select tolerances in relation to scaling, units, and required accuracy. Finite input validation does not prevent overflow during arithmetic on extreme magnitudes; the reference implementation targets moderate, sensibly scaled inputs and does not estimate a condition number.

from math import isfinite


def solve(a: list[list[float]], b: list[float], tol: float = 1e-12) -> list[float]:
    n = len(a)
    if n == 0 or len(b) != n or any(len(row) != n for row in a):
        raise ValueError("require a nonempty square system")
    if not isfinite(tol) or not 0 < tol < 1:
        raise ValueError("require 0 < tol < 1")
    if any(not isfinite(v) for row in a for v in row) or any(not isfinite(v) for v in b):
        raise ValueError("inputs must be finite")
    u, rhs = [row[:] for row in a], b[:]
    scale = max(abs(v) for row in a for v in row)
    for k in range(n):
        pivot = max(range(k, n), key=lambda i: abs(u[i][k]))
        if abs(u[pivot][k]) <= tol * scale:
            raise ValueError("singular or numerically rejected system")
        u[k], u[pivot] = u[pivot], u[k]
        rhs[k], rhs[pivot] = rhs[pivot], rhs[k]
        for i in range(k + 1, n):
            factor = u[i][k] / u[k][k]
            u[i][k] = 0.0
            for j in range(k + 1, n):
                u[i][j] -= factor * u[k][j]
            rhs[i] -= factor * rhs[k]
    x = [0.0] * n
    for i in range(n - 1, -1, -1):
        value = rhs[i]
        for j in range(i + 1, n):
            value -= u[i][j] * x[j]
        x[i] = value / u[i][i]
    return x


def sparse_matvec(columns: list[list[int]], values: list[list[float]],
                  x: list[float]) -> list[float]:
    if len(columns) != len(values):
        raise ValueError("row counts must match")
    out = [0.0] * len(columns)
    for i in range(len(columns)):
        if len(columns[i]) != len(values[i]):
            raise ValueError("row lengths must match")
        for column, value in zip(columns[i], values[i]):
            if not 0 <= column < len(x):
                raise ValueError("column index out of range")
            out[i] += value * x[column]
    return out


def residual_inf(a: list[list[float]], x: list[float], b: list[float]) -> float:
    result = 0.0
    for i, row in enumerate(a):
        value = 0.0
        for j, entry in enumerate(row):
            value += entry * x[j]
        result = max(result, abs(value - b[i]))
    return result


def show(values: list[float]) -> str:
    return "[" + ", ".join(f"{value:.6f}" for value in values) + "]"


def main() -> None:
    a = [[0.0, 2.0, 1.0], [1.0, 3.0, 0.0], [4.0, 0.0, 5.0]]
    b = [7.0, 7.0, 19.0]
    columns = [[1, 2], [0, 1], [0, 2]]
    values = [[2.0, 1.0], [1.0, 3.0], [4.0, 5.0]]
    x = solve(a, b)
    product = sparse_matvec(columns, values, x)
    error = residual_inf(a, x, b)
    print(f"nnz = {sum(len(row) for row in values)}")
    print(f"x = {show(x)}")
    print(f"sparse Ax = {show(product)}")
    print(f"residual_inf = {error:.3e}")
    print("residual check = " + ("PASS" if error <= 1e-10 else "FAIL"))


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <vector>

using Matrix = std::vector<std::vector<double>>;
using Vector = std::vector<double>;

Vector solve(const Matrix& a, const Vector& b, double tol = 1e-12) {
    const std::size_t n = a.size();
    if (n == 0 || b.size() != n) throw std::invalid_argument("invalid system shape");
    if (!std::isfinite(tol) || tol <= 0 || tol >= 1) throw std::invalid_argument("invalid tolerance");
    double scale = 0;
    for (const auto& row : a) {
        if (row.size() != n) throw std::invalid_argument("matrix must be square");
        for (double value : row) {
            if (!std::isfinite(value)) throw std::invalid_argument("inputs must be finite");
            scale = std::max(scale, std::abs(value));
        }
    }
    for (double value : b) if (!std::isfinite(value)) throw std::invalid_argument("inputs must be finite");
    Matrix u = a;
    Vector rhs = b;
    for (std::size_t k = 0; k < n; ++k) {
        std::size_t pivot = k;
        for (std::size_t i = k + 1; i < n; ++i) {
            if (std::abs(u[i][k]) > std::abs(u[pivot][k])) pivot = i;
        }
        if (std::abs(u[pivot][k]) <= tol * scale) throw std::runtime_error("singular or numerically rejected system");
        std::swap(u[k], u[pivot]);
        std::swap(rhs[k], rhs[pivot]);
        for (std::size_t i = k + 1; i < n; ++i) {
            double factor = u[i][k] / u[k][k];
            u[i][k] = 0;
            for (std::size_t j = k + 1; j < n; ++j) u[i][j] -= factor * u[k][j];
            rhs[i] -= factor * rhs[k];
        }
    }
    Vector x(n, 0);
    for (std::size_t index = n; index > 0; --index) {
        const std::size_t i = index - 1;
        double value = rhs[i];
        for (std::size_t j = i + 1; j < n; ++j) value -= u[i][j] * x[j];
        x[i] = value / u[i][i];
    }
    return x;
}

Vector sparse_matvec(const std::vector<std::vector<int>>& columns,
                     const Matrix& values, const Vector& x) {
    if (columns.size() != values.size()) throw std::invalid_argument("row counts must match");
    Vector out(columns.size(), 0);
    for (std::size_t i = 0; i < columns.size(); ++i) {
        if (columns[i].size() != values[i].size()) throw std::invalid_argument("row lengths must match");
        for (std::size_t k = 0; k < columns[i].size(); ++k) {
            int column = columns[i][k];
            if (column < 0 || static_cast<std::size_t>(column) >= x.size()) throw std::invalid_argument("column index out of range");
            out[i] += values[i][k] * x[static_cast<std::size_t>(column)];
        }
    }
    return out;
}

double residual_inf(const Matrix& a, const Vector& x, const Vector& b) {
    double result = 0;
    for (std::size_t i = 0; i < a.size(); ++i) {
        double value = 0;
        for (std::size_t j = 0; j < a[i].size(); ++j) value += a[i][j] * x[j];
        result = std::max(result, std::abs(value - b[i]));
    }
    return result;
}

void show(const char* label, const Vector& values) {
    std::cout << label << " = [" << std::fixed << std::setprecision(6);
    for (std::size_t i = 0; i < values.size(); ++i) {
        if (i != 0) std::cout << ", ";
        std::cout << values[i];
    }
    std::cout << "]\n";
}

int main() {
    const Matrix a = { {0, 2, 1}, {1, 3, 0}, {4, 0, 5}};
    const Vector b = {7, 7, 19};
    const std::vector<std::vector<int>> columns = { {1, 2}, {0, 1}, {0, 2}};
    const Matrix values = { {2, 1}, {1, 3}, {4, 5}};
    const Vector x = solve(a, b);
    const Vector product = sparse_matvec(columns, values, x);
    const double error = residual_inf(a, x, b);
    std::cout << "nnz = 6\n";
    show("x", x);
    show("sparse Ax", product);
    std::cout << "residual_inf = " << std::scientific << std::setprecision(3) << error << '\n';
    std::cout << "residual check = " << (error <= 1e-10 ? "PASS" : "FAIL") << '\n';
    return 0;
}

import java.util.Locale;

public class Main {
    public static double[] solve(double[][] a, double[] b, double tol) {
        int n = a.length;
        if (n == 0 || b.length != n) throw new IllegalArgumentException("invalid system shape");
        if (!Double.isFinite(tol) || tol <= 0 || tol >= 1) throw new IllegalArgumentException("invalid tolerance");
        double scale = 0;
        double[][] u = new double[n][];
        for (int i = 0; i < n; ++i) {
            if (a[i].length != n) throw new IllegalArgumentException("matrix must be square");
            u[i] = a[i].clone();
            for (double value : a[i]) {
                if (!Double.isFinite(value)) throw new IllegalArgumentException("inputs must be finite");
                scale = Math.max(scale, Math.abs(value));
            }
        }
        for (double value : b) if (!Double.isFinite(value)) throw new IllegalArgumentException("inputs must be finite");
        double[] rhs = b.clone();
        for (int k = 0; k < n; ++k) {
            int pivot = k;
            for (int i = k + 1; i < n; ++i) {
                if (Math.abs(u[i][k]) > Math.abs(u[pivot][k])) pivot = i;
            }
            if (Math.abs(u[pivot][k]) <= tol * scale) throw new IllegalArgumentException("singular or numerically rejected system");
            double[] row = u[k]; u[k] = u[pivot]; u[pivot] = row;
            double entry = rhs[k]; rhs[k] = rhs[pivot]; rhs[pivot] = entry;
            for (int i = k + 1; i < n; ++i) {
                double factor = u[i][k] / u[k][k];
                u[i][k] = 0;
                for (int j = k + 1; j < n; ++j) u[i][j] -= factor * u[k][j];
                rhs[i] -= factor * rhs[k];
            }
        }
        double[] x = new double[n];
        for (int i = n - 1; i >= 0; --i) {
            double value = rhs[i];
            for (int j = i + 1; j < n; ++j) value -= u[i][j] * x[j];
            x[i] = value / u[i][i];
        }
        return x;
    }

    public static double[] sparseMatvec(int[][] columns, double[][] values, double[] x) {
        if (columns.length != values.length) throw new IllegalArgumentException("row counts must match");
        double[] out = new double[columns.length];
        for (int i = 0; i < columns.length; ++i) {
            if (columns[i].length != values[i].length) throw new IllegalArgumentException("row lengths must match");
            for (int k = 0; k < columns[i].length; ++k) {
                int column = columns[i][k];
                if (column < 0 || column >= x.length) throw new IllegalArgumentException("column index out of range");
                out[i] += values[i][k] * x[column];
            }
        }
        return out;
    }

    public static double residualInf(double[][] a, double[] x, double[] b) {
        double result = 0;
        for (int i = 0; i < a.length; ++i) {
            double value = 0;
            for (int j = 0; j < a[i].length; ++j) value += a[i][j] * x[j];
            result = Math.max(result, Math.abs(value - b[i]));
        }
        return result;
    }

    private static void show(String label, double[] values) {
        System.out.print(label + " = [");
        for (int i = 0; i < values.length; ++i) {
            if (i != 0) System.out.print(", ");
            System.out.printf(Locale.ROOT, "%.6f", values[i]);
        }
        System.out.println("]");
    }

    public static void main(String[] args) {
        double[][] a = { {0, 2, 1}, {1, 3, 0}, {4, 0, 5}};
        double[] b = {7, 7, 19};
        int[][] columns = { {1, 2}, {0, 1}, {0, 2}};
        double[][] values = { {2, 1}, {1, 3}, {4, 5}};
        double[] x = solve(a, b, 1e-12);
        double[] product = sparseMatvec(columns, values, x);
        double error = residualInf(a, x, b);
        System.out.println("nnz = 6");
        show("x", x);
        show("sparse Ax", product);
        System.out.printf(Locale.ROOT, "residual_inf = %.3e%n", error);
        System.out.println("residual check = " + (error <= 1e-10 ? "PASS" : "FAIL"));
    }
}

Complexity, Correctness, and Invariants

For an $m\times n$ sparse row-list matrix with $s$ stored entries, multiplication takes $\Theta(m+s)$ time and $\Theta(m)$ space for the result, beyond the stored matrix and input vector. The matrix representation uses $\Theta(m+s)$ storage. Explicit zeros or duplicates count as stored entries for this runtime bound, even if a canonical mathematical nonzero count would be smaller. Dense multiplication takes $\Theta(mn)$ time.

For a successfully solved $n\times n$ system, the shown elimination code takes $\Theta(n^3)$ time and $\Theta(n^2)$ copied workspace. Pivot searches and back substitution each add $O(n^2)$ work. Its loops do not skip zero multipliers, so even a diagonal input follows the cubic update structure. Among accepted nonsingular systems of the same dimension, best, average, and worst asymptotic times are all cubic. Rejection in the first column can finish in $\Theta(n^2)$ time after validation and copying; an average including failures requires an input distribution.

In exact arithmetic, row exchanges and subtracting a multiple of one equation from another preserve the solution set. Each accepted pivot permits elimination below it. Earlier pivot columns stay zero below the diagonal because the pivot row is already zero in those columns. The final system is upper triangular and equivalent to the original, and back substitution uniquely determines its solution when every pivot is nonzero.

Floating-point arithmetic introduces rounding, so assess the residual and the conditioning separately. A useful scaled residual is

$$ \eta_\infty=\frac{\lVert Ax-b\rVert_\infty}{\lVert A\rVert_\infty\lVert x\rVert_\infty+\lVert b\rVert_\infty}, $$

where $\lVert A\rVert_\infty$ is the maximum absolute row sum. When the denominator is zero and the residual is zero, define this diagnostic as zero. A small scaled residual indicates a small normwise backward error; an ill-conditioned matrix can still amplify it into substantial solution error. The sample prints the unscaled residual for direct inspection.

Variants, Pitfalls, and Optimizations

Fill-in changes storage needs. The working entry at $(1,2)$ becomes $-5/4$ in the example even though that position was initially zero. A true sparse factorization must allocate new entries or predict their locations during symbolic analysis. Dropping fill-in changes the system unless the algorithm intentionally constructs an approximate factorization or preconditioner.

Separate ordering from pivoting. A fill-reducing ordering aims to make the factors sparse. Numerical pivoting chooses stable pivots during factorization. They serve different purposes and may interact. General sparse LU needs a specialized implementation; the original nonzero count alone does not determine its runtime or memory use.

Choose a solver for the matrix. Symmetric positive-definite systems can use Cholesky or conjugate gradients. General nonsymmetric systems may use sparse LU or methods such as GMRES with a suitable preconditioner. Iterative methods rely heavily on sparse matrix–vector products, but convergence depends on the matrix and preconditioning. The plain nested-array solver here is a reference for understanding elimination, not a scalable sparse factorization package.

Validate storage conventions. COO may contain repeated coordinates; define whether to sum them during assembly. Canonical CSR normally coalesces duplicates and often sorts column indices. Repeated row pointers represent empty rows, and the final pointer equals the stored-entry count. Converting a dense array by scanning it still costs $\Theta(mn)$; sparse multiplication's speed does not remove that construction cost.

Avoid explicit inversion for a solve. Solve using a factorization rather than forming $A^{-1}$. For many right-hand sides, retain the row permutation and LU factors. The teaching solver discards its factors, so calling it repeatedly repeats the cubic work. Iterative refinement can reuse a factorization while checking residuals against the original matrix.

Scale and diagnose thoughtfully. Swapping matrix rows without swapping the right-hand side corrupts the equations. Testing pivots against exact zero ignores near-singular cases, while a global threshold can reject badly scaled but invertible systems. Scaling rows or columns may help; rank-revealing QR or SVD is appropriate when rank or least-squares behavior is the actual question. A tolerance rejection alone does not distinguish inconsistent systems from systems with infinitely many solutions.

Practice Problems and Real-World Use

Choose storage and algorithms from the shape, pattern, and numerical properties of the matrix. The following exercises connect the worked operations to practical systems while providing small checks that can be verified by hand.

For verification, generate a small matrix and a known vector, compute the right-hand side independently, and recover the vector with the solver. Test row permutations, empty sparse rows, duplicate stored coordinates, malformed shapes, and nonfinite solve inputs. Include both a well-conditioned system and one that demonstrates why numerical tolerance and conditioning require separate interpretation.