Sparse Pseudoinverses

• Sparse pseudoinverses are generalized matrix inverses obtained by minimizing the ℓ1 norm to enforce optimal sparsity while preserving core Penrose properties.
• They are computed using convex optimization methods like ADMM and Douglas–Rachford splitting, enabling efficient least-squares recovery and scalability.
• Compared to the dense Moore–Penrose pseudoinverse, sparse pseudoinverses reduce storage and arithmetic costs with only a modest trade-off in numerical stability.

A sparse pseudoinverse is a generalized inverse of a matrix constructed to minimize an entrywise norm—predominantly the $\ell_1$ norm—as a surrogate for minimizing the number of nonzero entries, consistent with the fundamental inverse identity or further Penrose properties. This notion emerges as an alternative to the Moore–Penrose pseudoinverse, which is uniquely characterized by four Penrose equations but is typically dense and costly to apply, even for sparse or structured input matrices. Sparse pseudoinverses maintain algebraic well-posedness for inverse problems—such as least-squares—in high dimensions, while providing substantial gains in storage, arithmetic complexity, and interpretability, at modest cost to numerical stability and approximation quality (Dokmanić et al., 2017).

1. Definitions and Fundamental Properties

Let $A \in \mathbb{R}^{m \times n}$, typically with full column rank ($m < n$). A matrix $X\in\mathbb{R}^{n\times m}$ is a generalized inverse if $A X A = A$ (Penrose P1). For $1 \leq p < \infty$, the $\ell^p$-minimal inverse—termed the sparse pseudoinverse for $p=1$—is

$$\mathrm{SPI}_p(A) = \operatorname*{arg\,min}_{X}\|X\|_p \quad \text{s.t.} \quad A X A = A,$$

where $$\|X\|_p = \|\operatorname{vec}(X)\|_p = \left(\sum_{i,j}|X_{ij}|^p\right)^{1/p}$$. The sparse pseudoinverse is thus $A \in \mathbb{R}^{m \times n}$0, which solves

$A \in \mathbb{R}^{m \times n}$1

Existence follows from convexity of the objective and the affine constraint set. For generic $A \in \mathbb{R}^{m \times n}$2 (randomly drawn or in general position), the solution is unique and achieves optimal sparsity: each column of $A \in \mathbb{R}^{m \times n}$3 has exactly $A \in \mathbb{R}^{m \times n}$4 nonzeros, yielding $A \in \mathbb{R}^{m \times n}$5 nonzeros in total (Dokmanić et al., 2017).

A variety of structural constraints (e.g., reflexivity, symmetry, ah-symmetry) and block-support restrictions can be imposed to enforce additional Penrose properties (P2–P4) or target specific applications such as least-squares recovery, structured regression, or PDE-based approximation (Ponte et al., 2023, Ponte et al., 2024).

2. Comparison with the Moore–Penrose Pseudoinverse

The Moore–Penrose pseudoinverse ($A \in \mathbb{R}^{m \times n}$6) is uniquely characterized by the full Penrose system (P1–P4). However, $A \in \mathbb{R}^{m \times n}$7 is generically dense, requiring $A \in \mathbb{R}^{m \times n}$8 storage and $A \in \mathbb{R}^{m \times n}$9 flops to apply as a matrix–vector product. The sparse pseudoinverse, by contrast, requires only $m < n$0 nonzeros (for generic $m < n$1), yielding $m < n$2 application cost, which is a substantial improvement for $m < n$3.

For full-rank $m < n$4 with $m < n$5, the sparse pseudoinverse achieves the “sparsest possible pattern among all generalized inverses” and is generically unique (Dokmanić et al., 2017). In terms of numerical stability, sparse pseudoinverses increase the Frobenius norm (a conditioning measure) only marginally over $m < n$6; the mean-squared-error (MSE) penalty in noisy least-squares is controlled by a multiplicative factor $m < n$7 with $m < n$8 for random Gaussian $m < n$9, nearly matching that of $X\in\mathbb{R}^{n\times m}$0 over a broad range of $X\in\mathbb{R}^{n\times m}$1 (Dokmanić et al., 2017).

Inverse	Nonzeros	Apply Cost	Frobenius Norm ($X\in\mathbb{R}^{n\times m}$2)
$X\in\mathbb{R}^{n\times m}$3	$X\in\mathbb{R}^{n\times m}$4	$X\in\mathbb{R}^{n\times m}$5	$X\in\mathbb{R}^{n\times m}$6
$X\in\mathbb{R}^{n\times m}$7	$X\in\mathbb{R}^{n\times m}$8	$X\in\mathbb{R}^{n\times m}$9	$A X A = A$0 (slightly higher)

The trade-off is thus: sparse pseudoinverses are much sparser and faster, at the price of only a modest increase in $A X A = A$1 and resultant MSE (Dokmanić et al., 2017).

3. Norm-Minimization and Optimization Formulations

Sparse pseudoinverses are typically constructed by convex optimization, predominantly via:

• Entrywise $A X A = A$2-minimization (standard):

$A X A = A$3

This directly models entrywise sparsity.

• Mixed $A X A = A$4-norm minimization (row-sparsity):

$A X A = A$5

Minimized for row-sparse inverses—suited for applications where controllable row count is critical (Ponte et al., 2023, Ponte et al., 2024).

• Partial/Full Penrose constraints: Additional Penrose conditions (e.g., reflexivity, ah-symmetry) are incorporated to ensure desired properties:

$A X A = A$6

Depending on which are imposed, one obtains partial, reflexive, ah-symmetric, or symmetric sparse pseudoinverses (Fuentes et al., 2016, Fampa et al., 2020, Machado et al., 25 May 2026).

• SDP relaxations: To enforce quadratic Penrose constraints (e.g., $A X A = A$7), SDP "liftings" can be employed, but are limited to small-scale problems due to high computational complexity (Fuentes et al., 2016).

Efficient large-scale solvers for these convex programs use first-order methods such as ADMM and Douglas–Rachford splitting, exploiting block structure and efficient projection operators (Ponte et al., 2024, Machado et al., 25 May 2026).

4. Structural Results: Uniqueness, Sparsity, and Stability

Uniqueness and Support

For generic $A X A = A$8 (columns in general position and $A X A = A$9 not in low-dimensional spans), the $1 \leq p < \infty$0 minimization

$1 \leq p < \infty$1

admits unique solutions for every $1 \leq p < \infty$2; thus, $1 \leq p < \infty$3 is uniquely defined for almost all $1 \leq p < \infty$4 (Dokmanić et al., 2017).

Sparsity Bounds

For full column rank $1 \leq p < \infty$5, $1 \leq p < \infty$6 achieves optimal sparsity: each column contains exactly $1 \leq p < \infty$7 nonzeros, i.e., $1 \leq p < \infty$8 nonzeros in total, which is the minimal possible nonzero pattern for any generalized inverse (Dokmanić et al., 2017, Ponte et al., 2023).

Stability and Conditioning

The Frobenius norm of $1 \leq p < \infty$9 (conditioning) is tightly concentrated. For random Gaussian $\ell^p$0, the deviation of $\ell^p$1 from a deterministic $\ell^p$2 is exponentially small in $\ell^p$3. For $\ell^p$4 and $\ell^p$5, explicit formulas and concentration bounds are established (Dokmanić et al., 2017). The penalty incurred, compared to $\ell^p$6, in least-squares or denoising regimes is thus provably small.

5. Algorithmic Strategies and Computational Considerations

Convex Programming and First-Order Methods

Large-scale sparse pseudoinverse computation leverages splitting methods—primarily ADMM and Douglas–Rachford splitting—which allow efficient proximal iterations:

• The $\ell^p$7-minimization is handled via entrywise soft-thresholding.
• The affine constraint projection exploits the block structure derived via SVD or matrix symmetries for fast closed-form projection (Ponte et al., 2024, Machado et al., 25 May 2026).

Local-Search Heuristics

When strict row-sparsity or rank-minimality is required, combinatorial local-search on block-structured supports efficiently constructs ah-symmetric reflexive pseudoinverses: select a full-rank $\ell^p$8 submatrix, and set the inverse nonzero only in the corresponding rows/columns. Swapping operations are guided by local improvements in determinant or norm, yielding a $\ell^p$9- (or $p=1$0-) approximation to the true minimum $p=1$1 or $p=1$2 solution (Fampa et al., 2020, Ponte et al., 2023, Ponte et al., 2024).

Scalability and Practical Limits

State-of-the-art first-order algorithms (ADMM for $p=1$3 or $p=1$4 problems, DRS for symmetric inverses) scale to $p=1$5 in the thousands. General-purpose LP or SDP solvers are limited to $p=1$6 (LP) or $p=1$7 (SDP). Local-search methods are extremely fast for moderate ranks and yield strict support patterns suited to massive data or model compression (Ponte et al., 2024, Ponte et al., 2023, Machado et al., 25 May 2026).

6. Structured and Symmetric Sparse Pseudoinverses

Sparse pseudoinverse theory admits extensions to structured cases:

• Symmetric matrices: For symmetric $p=1$8, one seeks $p=1$9 with $$\mathrm{SPI}_p(A) = \operatorname*{arg\,min}_{X}\|X\|_p \quad \text{s.t.} \quad A X A = A,$$0 and symmetries in $$\mathrm{SPI}_p(A) = \operatorname*{arg\,min}_{X}\|X\|_p \quad \text{s.t.} \quad A X A = A,$$1. A new compact affine characterization $$\mathrm{SPI}_p(A) = \operatorname*{arg\,min}_{X}\|X\|_p \quad \text{s.t.} \quad A X A = A,$$2 enables efficient projection and DRS updates, maintaining symmetry and affine feasibility (Machado et al., 25 May 2026).
• Row- and block sparsity: Mixed norm minimization ($$\mathrm{SPI}_p(A) = \operatorname*{arg\,min}_{X}\|X\|_p \quad \text{s.t.} \quad A X A = A,$$3, $$\mathrm{SPI}_p(A) = \operatorname*{arg\,min}_{X}\|X\|_p \quad \text{s.t.} \quad A X A = A,$$4) directly controls row or column sparsity. For applications in model selection or interpretable regression, strict row-count can be imposed via combinatorial or nonconvex ADMM variants (Ponte et al., 2024, Ponte et al., 2023).

7. Domain Applications and Illustrative Examples

Sparse pseudoinverses are central in large-scale least-squares and regression with many right-hand sides, compressed sensing, PDE-based signal processing, statistics, and data science:

• In PDE-based sparse approximation (difference matrix case), the columns of the pseudoinverse form a dictionary for sparse signal representation. Orthogonal Matching Pursuit (OMP) and related greedy algorithms are then used to select a sparse subset of atoms (not pseudoinverse entries) for function approximation (Plonka et al., 2015).
• In operator-theoretic settings (reproducing kernel Hilbert spaces), the sparse pseudoinverse admits a sparse expansion as a greedy sum over RKHS kernels using POAFD or matching-pursuit strategies, guaranteeing $$\mathrm{SPI}_p(A) = \operatorname*{arg\,min}_{X}\|X\|_p \quad \text{s.t.} \quad A X A = A,$$5 convergence (Qian, 2019).
• Large-scale machine learning applications benefit from sparse or low-rank approximations (e.g., FastPI–incremental SVD for sparse matrices), ensuring scalability with minimal accuracy loss (Jung et al., 2020).

• "Beyond Moore-Penrose Part II: The Sparse Pseudoinverse" (Dokmanić et al., 2017)
• "Good and Fast Row-Sparse ah-Symmetric Reflexive Generalized Inverses" (Ponte et al., 2024)
• "Sparse pseudoinverses via LP and SDP relaxations of Moore-Penrose" (Fuentes et al., 2016)
• "Trading off 1-norm and sparsity against rank for linear models using mathematical optimization: 1-norm minimizing partially reflexive ah-symmetric generalized inverses" (Fampa et al., 2020)
• "On computing sparse generalized inverses" (Ponte et al., 2023)
• "Sparse symmetric generalized inverses for sparse symmetric matrices" (Machado et al., 25 May 2026)
• "Experimental analysis of local searches for sparse reflexive generalized inverses" (Fampa et al., 2020)
• "Pseudo-inverses of difference matrices and their application to sparse signal approximation" (Plonka et al., 2015)
• "Sparse Solutions for Inverse Problems in Reproducing Kernel Hilbert Spaces" (Qian, 2019)
• "Fast and Accurate Pseudoinverse with Sparse Matrix Reordering and Incremental Approach" (Jung et al., 2020)