Factorization-Free Matrix Computation

• Factorization-free matrix computation is a set of techniques that bypass explicit matrix factorizations to maintain exact arithmetic and streamline high-dimensional computations.
• It employs methods such as Gauss–Bareiss reduction, matrix-free operator evaluations, and tensor product factorizations to control growth and enhance numerical stability.
• The topic also explores pivoting strategies, localized divide-and-conquer algorithms, and implicit regularization to achieve efficient, scalable, and robust solutions.

Factorization-free matrix computation encompasses a spectrum of exact and approximate linear algebraic methods that deliberately avoid traditional matrix factorization steps—such as LU, QR, or explicit spectral decompositions—whether for reasons of computational efficiency, numerical stability, exact arithmetic, or scalability. These approaches circumvent or replace direct matrix factorization by strategies including fraction-free arithmetic, localized divide-and-conquer algorithms, implicit regularization, matrix-free operator applications, and exploitation of tensor product structure for modular high-dimensional computations.

1. Fraction-Free Decomposition via Gauss–Bareiss Reduction

The Gauss–Bareiss reduction algorithm provides a central paradigm for fraction-free computation. Contrary to classical Gaussian elimination over fields (which can rapidly introduce fractions even for integral input), Bareiss reduction maintains all entries within the original domain by exact division at each step. Given a matrix $A$, this leads to an LU-like decomposition:

$A = P_w \cdot L \cdot D^{-1} \cdot U \cdot P_c$

where $P_w$ and $P_c$ are permutation matrices, $L$ and $U$ are lower and upper triangular, and $D$ is a diagonal matrix capturing the products of pivotal divisors. By organizing arithmetic to cancel emerging denominators as soon as possible, the algorithm precludes fraction growth and maintains exact arithmetic throughout.

Fraction-free QR algorithms, notably those based on Gram–Schmidt processes, yield a reduced QR-style decomposition without introducing fractions or irrationalities. The reduced factor $\Theta$ satisfies $\Theta^T \Theta = D$ for a diagonal $D$. Notably, for square, full-rank matrices, the last column of $\Theta$ manifests a systematic factor proportional to $\det A$, as shown by explicit determinantal formulae.

2. Systematic and Statistical Common Factors in Decomposition

In fraction-free LU and QR decompositions, two classes of common factors arise in the output:

• Systematic factors: Determined by the elimination algorithm structure. QR produces a $\det A$ factor in the $\Theta$ matrix's final column, while LU decompositions naturally accumulate pivot products as consistent row or column factors. These align closely with the determinantal divisors appearing in Smith normal forms. Specifically, for $A = L D^{-1} U$, the product $d^*_k = \prod_{j=1}^{k} d_j$ divides every entry of row $k$ in $U$ and column $k$ in $L$.
• Statistical factors: Emergent from specific input data properties. For integer matrices, elimination steps can induce new accidental common divisors through cross-multiplications. The probability of statistical factor occurrence in an elimination step is approximately 27% for pairs of integers, with empirical observations indicating that 40% of predicted row factors actually manifest in practice.

This taxonomy situates systematic factors as traceable to algebraic invariants, while statistical factors are data-induced and estimated probabilistically.

3. Pivoting Strategies and Output Size Control

Fraction-free algorithms eschew pivoting for numerical stability but employ it to control intermediate and output size. Several pivoting strategies are analyzed:

• Largest: Maximizes the pivot choice by absolute value or degree; this produces more voluminous intermediate entries.
• Smallest: Minimizes the pivot via the same criteria, consistently yielding more compact outputs (fewer digits or polynomial coefficients) and often reducing computational complexity.
• First: Selects the earliest nonzero, offering simplicity but suboptimal control of growth.
• Factors: Chooses pivots by minimal number of prime divisors, which theoretically minimizes output but is computationally prohibitive.

Experimental data confirm that "smallest" pivoting strategies lead to more efficient decompositions for integer and polynomial matrices, reducing total digit count and the frequency of row GCD factors in $U$ compared to "largest" or "first" strategies.

4. Localized Divide-and-Conquer Methods for Inverse Factorization

For Hermitian positive definite matrices with exponential decay structure (e.g., those arising in quantum chemistry), recursive localized inverse factorization is a factorization-free paradigm (Rubensson et al., 2018). The method partitions $S$ into submatrices (e.g., $A$ and $C$), solves $Z_A$ and $Z_C$ such that $Z_A^* A Z_A = I$ and analogously for $C$, and uses a block-diagonal initial guess $Z_0$. Iterative refinement "glues" subproblem solutions by correcting only near the interface ("cut"), reducing computation to $O($cut size$)$ under locality assumptions. Iterative updates follow polynomial iterations akin to Newton–Schulz methods, with convergence rates dependent on the condition number of $S$.

The method is fully parallelizable, with subproblems solved independently and combined with negligible overhead for large matrices where the cut size grows sublinearly. Theoretical results provide exponential decay bounds for corrections away from the cut, and empirical demonstrations include lattice adjacency matrices and basis set overlaps from Hartree–Fock and Kohn–Sham codes.

5. Matrix-Free Operator Evaluation and Jacobian Chaining

Efficient large-scale computations often require matrix-free approaches—in which operators such as Jacobians are never explicitly assembled or factorized, but rather applied directly to vectors (Naumann, 11 Apr 2024). For differentiable programs comprising sequences of subprograms, the Jacobian chain product

$F' = F^{(q)'} \cdot F^{(q-1)'} \cdots F^{(1)'}$

is computed via chained tangent (forward mode) and adjoint (reverse mode) applications, avoiding explicit matrix construction. The matrix-free protocol is formulated dynamically, balancing fused multiply–add cost and tape memory with checkpointing to meet allocation constraints; bracketing strategies are selected via dynamic programming to minimize resources.

Numerical evidence shows reductions in computational work—sometimes more than an order of magnitude compared to classical approaches—validating the matrix-free method's efficiency, particularly under tight memory regimes for large simulation codes.

6. Matrix-Free Factorization Using Tensor Product Structure (Stabilization in CutFEM)

In the context of high-order finite element methods for problems with evolving geometries (CutFEM), matrix-free strategies leverage the tensor product structure of basis functions and operators (Wichrowski, 28 Feb 2025). Ghost penalty stabilization, crucial for small cut cells, is implemented without assembling the global matrix. The operator is decomposed into Kronecker products of 1D mass and penalty matrices: $\mathcal{G}_{F,k} = M^h \otimes M^h \otimes G^h_k$, where $G^h_k$ penalizes normal derivative jumps.

Sequential sum-factorization and precomputed 1D matrices reduce all multi-dimensional product evaluations to efficient 1D operations, yielding computational complexity $O(k^{d+1})$ for degree-$k$ elements in $d$ dimensions. Implementation within the deal.II library exploits batch vectorization and MPI parallelism. This matrix-free tensor product factorization is central for scalable, parallelizable CutFEM simulations.

7. Algorithmic Regularization in Model-Free Overparametrized Matrix Factorization

In unconstrained, nonconvex asymmetric matrix factorization tasks (with arbitrary overparametrization), global optima tend to overfit, interpolating noise. However, gradient descent initialized with sufficiently small random factors and subjected to early stopping implicitly regularizes the solution, sequentially recovering principal components without explicit regularization (Jiang et al., 2022).

For the objective $f(F, G) = \frac{1}{2} \| FG^T - X \|^2$, where $X$ is the observed matrix, iterates from small random initialization exhibit geometric growth only in directions of leading singular values, approaching the best rank-$r$ approximation $X_r$ when stopped appropriately. The iteration complexity to achieve error $\epsilon$ scales as $O(\log(1/\epsilon))$, nearly dimension-free except for logarithmic dependencies, with explicit bounds relating the stepsize, initialization, and singular value gaps. Empirical results confirm theoretical predictions of implicit regularization and highlight the importance of initialization and stopping strategy.

Summary

Factorization-free matrix computation comprises a set of coordinated strategies to circumvent traditional matrix factorization steps, maintaining exact arithmetic, exploiting locality and parallelism, leveraging algorithmic differentiation, and utilizing operator structure for high-dimensional problems. These methods—spanning fraction-free arithmetic, localized inverse factorizations, chain rule operator applications, and tensor product decomposition—advance the efficiency, scalability, and numerical robustness of matrix computations across exact, approximate, and simulation-based settings. Experimental, theoretical, and applied evidence across domains including symbolic computation, electronic structure calculations, large-scale numerical simulation, and high-order finite elements affirms the significance and versatility of the factorization-free approach.