PyTorch PR #200536: Why Small Matrices Returned Zero Pseudoinverses on MPS

For the matrix A = 1e-8 * I, the pseudoinverse is 1e8 * I. PyTorch PR #200536 reports that its MPS backend instead returned an all-zero matrix. Solving A X = I through torch.linalg.lstsq returned zeros as well, while the CPU results contained the expected 100000000 diagonal.

The failure came from PyTorch's Metal SVD kernel. It calculated nonzero singular values, then discarded their corresponding vector directions during normalization because the values were smaller than a fixed machine-epsilon threshold.

The patch changes that normalization condition from sigma > eps to sigma > 0.0f. It also adds tests that check orthonormality, reconstruction, pseudoinverses, and least-squares solutions at a small input scale.

The useful distinction is between the size of a matrix and whether its information has become numerically insignificant. This example has four equal singular values and condition number one. Scaling an identity matrix down did not make it rank-deficient.

The reported result was an arithmetic failure

The PR includes this reproducer:

import torch

for device in ("cpu", "mps"):
    eye = torch.eye(4, dtype=torch.float32, device=device)
    a = 1e-8 * eye
    print(device, "pinv:", torch.linalg.pinv(a).cpu())
    print(device, "lstsq:", torch.linalg.lstsq(a, eye).solution.cpu())

The expected answer follows directly from the input. Multiplying 1e-8 * I by 1e8 * I gives I, so the inverse, pseudoinverse, and solution against the identity agree.

An all-zero answer gives A X = 0. It cannot satisfy the requested system.

The MPS output in the PR is the contributor's reported result. I did not run this reproducer on Apple hardware. The locally executed checks described below use CPU PyTorch and an isolated model of the changed normalization rule.

Epsilon is not the smallest nonzero float

The kernel defines eps using the epsilon of float, approximately 1.1920929e-7.

That quantity describes floating-point spacing around one. It is not a universal boundary below which numbers must be treated as zero. A float32 can represent 1e-8 as a normal, nonzero value; even its square, near 1e-16, is well within the normal range.

The default pseudoinverse policy is also scale-aware. PyTorch documents a cutoff based on the largest singular value multiplied by a relative tolerance, together with any specified absolute tolerance. For a 4-by-4 float32 matrix with default settings, the relative tolerance is 4 * eps.

For this input, that gives a cutoff near:

1e-8 * 4 * 1.1920929e-7
= 4.7683716e-15

Every singular value, approximately 1e-8, is comfortably above it. The unexpected zeros originated before this legitimate rank-selection step could use the decomposition correctly.

The kernel kept the singular value and erased its vector

The reduced singular value decomposition expresses a matrix as:

A = U diag(S) Vh

Here Vh is the conjugate transpose of V. The columns of U and V must be orthonormal.

In the reviewed Jacobi kernel, the transformed columns contain both a direction and a magnitude. If a column is c = sigma * u, recovering its unit direction requires:

u = c / sigma

The kernel computes each column norm, sorts the resulting singular values, and writes sigma into S. It then computes the reciprocal used to emit the vector:

// Before PR #200536:
float inv = sigma > eps ? (1 / sigma) : 0.0f;

For sigma = 1e-8, the comparison fails. Multiplying the column by inv writes zeros into the vector output even though the singular-value output remains nonzero.

This creates an internally inconsistent decomposition: S says that a direction has nonzero magnitude, while its vector contains no direction at all. Inspecting singular values alone can miss the problem.

A later rank check did not repair these columns

The same kernel already has a separate step for completing the basis of numerically null directions. Its threshold is relative:

threshold = eps * largest_singular_value * working_row_count

It completes the trailing columns classified as numerically zero.

For the small scaled identity, all singular values remain above that relative threshold. The computed rank is full, so the completion step has no trailing null columns to replace. It leaves the incorrectly zeroed vectors untouched.

The two stages were therefore applying incompatible rules. Normalization rejected a direction because its magnitude was below an absolute epsilon, while basis completion recognized that same direction as significant relative to the matrix.

The patch removes the absolute rejection from normalization. It leaves the separate relative rank and basis-completion logic in place.

One comparison fixes both matrix orientations

The changed line is:

float inv = sigma > 0.0f ? (1 / sigma) : 0.0f;

At this output stage, a positive computed singular value receives its reciprocal. An exact zero still takes the zero branch, avoiding division by zero before the existing basis-completion logic runs.

The location covers two orientations. The MPS host implementation runs the Jacobi kernel with at least as many rows as columns. For a wide matrix, it processes the conjugate transpose and swaps the output roles.

For the ordinary orientation, the normalized column is written into U. For the transposed orientation, it is written into the corresponding Vh row, with conjugation.

That is why the regression uses both (8, 4) and (4, 8) matrices. The same predicate could erase a different returned factor depending on the shape.

Both higher-level operations depended on the broken factors

The pseudoinverse is constructed from the SVD factors:

A+ = V diag(S+) Uᴴ

Taking reciprocals of valid singular values cannot recover a vector that the SVD has already replaced with zero.

The MPS least-squares implementation exposes the same dependency directly. It obtains U, S, and Vh through linalg_svd, applies a relative singular-value cutoff, and computes:

tmp = (Uᴴ B) scaled by reciprocal singular values
solution = V tmp

A zeroed U or V factor contaminates the solution regardless of whether the later cutoff is correct. This explains why one kernel edit addresses the two public API symptoms in the reproducer.

The regression tests identities, not a particular basis

The new test uses rectangular identity matrices scaled by 1e-8, with float32 and complex64 storage. It requests a reduced SVD and checks five properties for each shape and dtype combination.

The first two are:

Uᴴ U = I
Vh Vhᴴ = I

These catch erased vectors immediately. The remaining checks verify rescaled reconstruction, the pseudoinverse, and the least-squares solution against an identity right-hand side.

It does not require the returned vectors to match a particular CPU basis element by element. Singular vectors are not unique, especially when singular values repeat. Orthonormality and reconstruction test the mathematical result without rejecting another valid basis.

Although the test includes complex storage, its identity inputs have zero imaginary parts. It exercises the complex kernel specialization but does not independently stress arbitrary complex phases.

Rescaling prevents an all-zero result from passing

A direct reconstruction check on a matrix whose entries are only 1e-8 can be misleading if it uses an absolute tolerance of 1e-5.

A completely zero reconstruction differs by only 1e-8, which is smaller than that tolerance. The test could accept a 100% relative error.

The added regression instead divides the singular values by the input scale before reconstruction:

((U * (S / scale).unsqueeze(-2)) @ Vh).cpu()

It compares that result with the unscaled identity. The expected diagonal is now one. A zero reconstruction differs by one and cannot hide beneath the tolerance.

It similarly multiplies the pseudoinverse and least-squares solution by scale before comparing them with the unscaled reference pseudoinverse.

This is a particularly useful testing pattern for small-scale numerical failures: bring the expected result back to a meaningful comparison scale instead of making the tolerance smaller by guesswork.

What was checked locally

On Linux with PyTorch 2.10.0+cpu, the square reference returned condition number one and the expected 100000000 diagonal for both operations.

I also applied only the old and new vector-emission predicates to correct CPU SVD factors. For the four shape/dtype combinations, the old predicate produced a maximum rescaled reconstruction error of one; the new predicate produced zero error. This models the disputed decision, not the complete Metal kernel.

The local tolerance check confirmed that an all-zero reconstruction passes an unscaled atol=1e-5 comparison and fails the rescaled comparison.

There was also a limitation worth retaining: the installed older CPU runtime's default lstsq failed the wide-matrix checks for both dtypes. Explicit CPU driver="gelsd" checks passed. The cause was not investigated, and neither result is treated as execution evidence for the MPS patch. The Research Report preserves those failures separately.

Scope of the repair

This change targets PyTorch's native Metal Jacobi SVD path for float32 and complex64. The reviewed host code falls back to CPU for matrices that exceed its threadgroup staging-memory limit, so an MPS tensor does not imply that every shape uses this kernel.

The added cases are small, full-rank, reduced-SVD examples at one scale. They do not establish accuracy across the entire floating-point range, genuinely rank-deficient matrices, batched layouts, or SVD gradients. No performance benchmark is provided.

The implementation change is consistent with the reported finite small-scale failure, but this review did not build the landed commit or execute the upstream MPS tests. A deployed build should be checked with orthonormality and scale-aware reconstruction, not just its package version or singular values.

The central rule is simple: a normalization guard must not erase a nonzero direction merely because the whole matrix is expressed in small units.