healpix_resample.psf.least_squares_cg

healpix_resample.psf.least_squares_cg(M, MT, y, x_ref, x0, max_iter=200, tol=1e-06, damp=0.0, verbose=True, weight=None)

Solve for delta in a damped least-squares problem without forming dense matrices.

M and MT are not Euclidean transposes of one another (each is normalized against a different axis of the raw weight matrix), but MT is exactly the adjoint of M with respect to a pair of weighted inner products: <.,.>_Dy on the sample space (weight = per-sample row-sum used to normalize MT) and <.,.>_Dx on the HEALPix-cell space (weight = per-cell column-sum used to normalize M, passed here as weight). Concretely MT = Dx @ M.T @ Dy^-1 (in this module’s row-vector convention). This solves the stationarity condition of the weighted least-squares problem

delta_hat = argmin_delta || delta @ MT - r_ref ||^2_Dy + damp * || delta ||^2_Dx

where r_ref = y - x_ref @ MT is the sample-space residual, which is exactly:

(MT-then-M + damp*I) delta = (y - x_ref @ MT) @ M

i.e. the same linear system as a naive (unweighted) Tikhonov normal equation, but its correct interpretation – and the correct inner product for the Conjugate Gradient solver below – uses Dx = weight (see conjugate_gradient()).

Shapes:

M : (N, K) sparse CSR MT : (K, N) sparse CSR y : (B, N) x_ref : (B, K) delta : (B, K) weight : (K,) or None – Dx, the per-cell weight columns of M were normalized by; None falls back to the (formally unjustified) Euclidean inner product.