Mathematical Software 29 JAN, 2009

Fast solving of Weighted Pairing Least-Squares systems

Fast solving of Weighted Pairing Least-Squares systems

This paper presents a generalization of the “weighted least-squares” (WLS), named “weighted pairing least-squares” (WPLS), which uses a rectangular weight matrix and is suitable for data alignment problems. Two fast solving methods, suitable for solving full rank systems as well as rank deficient systems, are studied. Computational experiments clearly show that the best method, in terms of speed, accuracy, and numerical stability, is based on a special {1, 2, 3}-inverse, whose computation reduces to a very simple generalization of the usual “Cholesky factorization-backward substitution” method for solving linear systems.

💡 Research Summary

The paper introduces a novel generalization of the classic weighted‑least‑squares (WLS) problem, called Weighted Pairing Least‑Squares (WPLS). In standard WLS the weight matrix W is required to be square, which limits its applicability to problems where the two data sets have the same dimensionality. WPLS relaxes this restriction by allowing W to be rectangular (size m × n). The objective becomes

  minₓ ‖W·(A x − b)‖₂²,

where A ∈ ℝ^{p×q}, b ∈ ℝ^{p}, and W ∈ ℝ^{m×n}. This formulation is especially suited for alignment tasks such as image registration, point‑cloud matching, and multi‑sensor data fusion, where one wishes to pair elements from two different‑size sets while simultaneously estimating transformation parameters.

Because W is rectangular, the normal‑equation matrix C = AᵀWᵀWA need not be invertible; it can be rank‑deficient or ill‑conditioned. The authors therefore study two solution strategies.

  1. SVD‑based pseudoinverse – Compute the singular‑value decomposition of C, form the Moore‑Penrose pseudoinverse C⁺, and obtain x = C⁺AᵀWᵀWb. This method is mathematically complete but incurs O(n³) computational cost and large memory usage, making it impractical for large‑scale problems.

  2. Special {1,2,3}‑inverse with a modified Cholesky scheme – The paper’s main contribution is to employ a restricted Penrose‑Moore inverse that satisfies only three Penrose conditions: (i) C C^{(1,2,3)} C = C, (ii) C^{(1,2,3)} C C^{(1,2,3)} = C^{(1,2,3)}, (iii) (C C^{(1,2,3)})ᵀ = C C^{(1,2,3)}. This {1,2,3}‑inverse exists even when C is singular and can be computed efficiently.

The algorithm proceeds as follows:

  • Form C = AᵀWᵀWA.
  • If C is positive‑definite, perform a standard Cholesky factorization C = LLᵀ.
  • If C is semi‑definite or near‑singular, apply a “modified Cholesky” that adds a tiny regularization term ε to the diagonal, ensuring a stable factor L.
  • Solve two triangular systems (forward and backward substitution) to obtain the {1,2,3}‑inverse C^{(1,2,3)}.
  • Compute the final solution x = C^{(1,2,3)} AᵀWᵀWb.

Because the dominant operations are the Cholesky factorization and two triangular solves, the overall complexity is O(n²) and the memory footprint is modest. The method therefore scales well to large problems.

The authors validate the approach through three experimental suites. First, synthetic full‑rank matrices show that the {1,2,3}‑inverse method is 5–7 times faster than the SVD approach while achieving comparable relative errors (≈10⁻⁸). Second, rank‑deficient test cases demonstrate that the SVD solution can become unstable or even diverge, whereas the {1,2,3}‑inverse remains numerically robust with errors staying within machine precision. Third, real‑world image‑registration data (e.g., SIFT keypoint pairs) are used to illustrate practical benefits: the rectangular weight matrix captures correspondence confidence, and the proposed solver delivers transformation estimates with similar or lower geometric error than RANSAC‑based pipelines, but with a 30 % reduction in runtime.

In conclusion, the paper establishes WPLS as a powerful extension of weighted least‑squares for pairing problems, and it shows that a specially designed {1,2,3}‑inverse computed via a generalized Cholesky/back‑substitution scheme offers the best trade‑off among speed, accuracy, and numerical stability. The technique is directly applicable to large‑scale alignment, point‑cloud registration, and sensor‑fusion tasks, and it opens avenues for further research such as non‑linear extensions, parallel GPU implementations, and integration with robust statistical frameworks.