Indefinite Least Squares Problem

• ILS is defined as minimizing the quadratic form (b-Ax)ᵀJ(b-Ax) with a signature matrix, which introduces negative weights and distinguishes it from ordinary least squares.
• Its block reformulations and splitting iterations enable the development of efficient preconditioners and iterative solvers, such as hyperbolic QR and randomized variants.
• Robust perturbation analysis and backward error estimates are crucial for understanding the stability and sensitivity of both constrained and unconstrained ILS problems.

Indefinite least squares (ILS) is the problem of minimizing an indefinite quadratic form of the residual,

$\min_{x\in\mathbb R^n}(b-Ax)^T J (b-Ax),$

where $A\in\mathbb R^{m\times n}$, $b\in\mathbb R^m$, and $J=\operatorname{diag}(I_p,-I_q)$ is a signature matrix with $p+q=m$. When $p=0$ or $q=0$, the problem reduces to the standard linear least squares problem; for a genuine ILS problem one assumes $A^T J A$ is positive definite, which guarantees a unique minimizer and yields the normal equations $A^T J A\,x=A^T J b$ (Diao et al., 2016). Relative to ordinary least squares, the formal change is the replacement of $A^T A$ by $A\in\mathbb R^{m\times n}$0, so some rows of $A\in\mathbb R^{m\times n}$1 enter with negative weight; this algebraic modification drives a distinct theory of block formulations, stationary iterations, Krylov preconditioning, and perturbation analysis (Grcar, 2010).

1. Algebraic formulation

The standard notation partitions the data conformally with the signature: $A\in\mathbb R^{m\times n}$2 with $A\in\mathbb R^{m\times n}$3, $A\in\mathbb R^{m\times n}$4, and $A\in\mathbb R^{m\times n}$5. Under this partition,

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

so the normal equations become

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

Several works assume $A\in\mathbb R^{m\times n}$8 has full column rank and $A\in\mathbb R^{m\times n}$9 is SPD; under these hypotheses the minimizer is unique (Zhang et al., 2022, Salkuyeh, 22 Jul 2025).

Expanding the objective,

$b\in\mathbb R^m$0

shows that the definiteness of $b\in\mathbb R^m$1 governs solvability. In the unconstrained case this is the basic well-posedness condition; in equality-constrained variants, positivity is instead required on $b\in\mathbb R^m$2 or $b\in\mathbb R^m$3, depending on the notation adopted in a given formulation (Li et al., 7 May 2026, Wang et al., 2016, Diao et al., 2018).

2. Equivalent block formulations

A recurring feature of ILS is its reformulation as a $b\in\mathbb R^m$4 block linear system by introducing residual variables. Let $b\in\mathbb R^m$5, partitioned as $b\in\mathbb R^m$6. One representation writes

$b\in\mathbb R^m$7

and this system is equivalent to the normal equations when $b\in\mathbb R^m$8 is SPD (Shaldehi et al., 28 Feb 2026).

Other papers use algebraically equivalent block systems that emphasize different elimination patterns or sparsity structures.

Source	Unknowns	Block coefficient matrix
(Shaldehi et al., 28 Feb 2026)	$b\in\mathbb R^m$9	$J=\operatorname{diag}(I_p,-I_q)$0
(Salkuyeh, 22 Jul 2025)	$J=\operatorname{diag}(I_p,-I_q)$1	$J=\operatorname{diag}(I_p,-I_q)$2
(Li et al., 23 May 2025)	$J=\operatorname{diag}(I_p,-I_q)$3	$J=\operatorname{diag}(I_p,-I_q)$4

In the parameterized formulation, $J=\operatorname{diag}(I_p,-I_q)$5 and $J=\operatorname{diag}(I_p,-I_q)$6, while in the sparse reformulation the middle diagonal block is $J=\operatorname{diag}(I_p,-I_q)$7, which becomes the target of a diagonal correction in later preconditioners (Salkuyeh, 22 Jul 2025, Li et al., 23 May 2025). This suggests that the block reformulation is not merely a change of variables: it is the structural basis for most modern ILS preconditioners and convergence proofs.

3. Direct solvers and perturbation questions

A classical direct approach is the hyperbolic QR factorization. In this setting one factors

$J=\operatorname{diag}(I_p,-I_q)$8

with $J=\operatorname{diag}(I_p,-I_q)$9 upper triangular and $p+q=m$0 $p+q=m$1-orthogonal, then solves the resulting triangular system after forming $p+q=m$2. Algorithmically this parallels Euclidean QR, but the orthogonal transformations are replaced by hyperbolic rotations or reflections that preserve the bilinear form induced by $p+q=m$3 (Grcar, 2010).

A second direct strategy is the QR–Cholesky method due to Chandrasekaran–Gu–Sayed. One first computes a thin hyperbolic QR factorization $p+q=m$4, then forms

$p+q=m$5

by Cholesky factorization, and finally solves

$p+q=m$6

Because the triangular factors $p+q=m$7 and $p+q=m$8 are already available, backward-error and condition estimates can be integrated into the same factor-solve framework at low extra cost (Diao et al., 2016).

The forward-stability theory of direct methods is not entirely straightforward. Grcar showed that a perturbation bound used by Bojanczyk, Higham, and Patel can produce arbitrarily large overestimates for all perturbations of some ILS problems, and therefore the hyperbolic QR factorization algorithm is not proved to be forward stable by that bound (Grcar, 2010). The underlying mechanism is cancellation in Jacobian-based expressions for the sensitivity with respect to $p+q=m$9, so a bound obtained through triangle inequalities may be much larger than the attainable first-order variation. This is one of the main technical cautions in the ILS literature.

4. Splitting iterations and randomized variants

A basic stationary iteration is the splitting (SP) method, obtained from

$p=0$0

by rewriting it as

$p=0$1

Since the eigenvalues of $p=0$2 satisfy $p=0$3, the SP iteration converges linearly for any initial vector. The same work develops two randomized variants, SP-RK-RGS and SP-SCD, by replacing the solve with $p=0$4 by randomized Kaczmarz, randomized Gauss–Seidel, or sampling coordinate descent steps; in the reported tests SP-SCD was the fastest, with speed-up $p=0$5 in one large random test and $p=0$6 in a Minkowski-space case relative to USSOR (Zhang et al., 2022).

A parameterized block splitting introduces

$p=0$7

with iteration matrix $p=0$8. Its convergence is characterized exactly: $p=0$9 and the optimal parameter is

$q=0$0

Away from $q=0$1, the eigenvalues of the preconditioned matrix $q=0$2 cluster in a circle centered at $q=0$3, which is favorable for GMRES (Salkuyeh, 22 Jul 2025).

A more recent development is the double-splitting (DS) iteration for the normal equations. With

$q=0$4

the two-step recurrence is

$q=0$5

By analyzing the companion matrix

$q=0$6

one obtains $q=0$7 for every $q=0$8, so the method converges unconditionally. In the reported dense and sparse experiments the DS method reached the target residual in two iterations and used less CPU time than SP, GSP, and ADI (Li et al., 7 May 2026).

5. Preconditioned GMRES and spectral clustering

The most active recent line of work concerns block-splitting preconditioners for GMRES. In the inexact block-splitting (IBS) framework, the exact middle block $q=0$9 is replaced by a well-conditioned SPD approximation $A^T J A$0. For the four splittings $A^T J A$1, the stationary iteration is

$A^T J A$2

Convergence conditions are given explicitly: $A^T J A$3 for IBS1 and IBS3, and

$A^T J A$4

for IBS2 and IBS4. Under these hypotheses $A^T J A$5, every eigenvalue $A^T J A$6 of $A^T J A$7 satisfies $A^T J A$8, $A^T J A$9 has large geometric multiplicity, and for IBS2 and IBS4 the non-unit eigenvalues lie on the positive real line in $A^T J A\,x=A^T J b$0. The minimal polynomial has degree at most $A^T J A\,x=A^T J b$1, so in exact arithmetic GMRES terminates in at most $A^T J A\,x=A^T J b$2 steps (Shaldehi et al., 28 Feb 2026).

A different accelerated preconditioner starts from the sparse $A^T J A\,x=A^T J b$3 system

$A^T J A\,x=A^T J b$4

and replaces the central zero block by $A^T J A\,x=A^T J b$5: $A^T J A\,x=A^T J b$6 If

$A^T J A\,x=A^T J b$7

then the fixed-point iteration under the splitting converges, the eigenvalues of $A^T J A\,x=A^T J b$8 are real, and as $A^T J A\,x=A^T J b$9 they cluster at $A^T A$0. In the reported GMRES experiments this preconditioner required $A^T A$1–$A^T A$2 iterations and the smallest CPU time across all tests (Li et al., 23 May 2025).

The numerical evidence in these papers is consistent with the spectral results. For the TOLS340 test in the IBS study, IBS2 and IBS4 required $A^T A$3 iterations and $A^T A$4 seconds, whereas BS2 and BUT required $A^T A$5 and $A^T A$6 iterations with $A^T A$7 and $A^T A$8 seconds, respectively; all tests there used FGMRES with inexact inner CG solves of tolerance $A^T A$9 and outer residual tolerance $A\in\mathbb R^{m\times n}$00 (Shaldehi et al., 28 Feb 2026). In the accelerated-preconditioner study, TOLS340 was solved in $A\in\mathbb R^{m\times n}$01 GMRES iterations and $A\in\mathbb R^{m\times n}$02 seconds with $A\in\mathbb R^{m\times n}$03, compared with $A\in\mathbb R^{m\times n}$04 iterations for BS$A\in\mathbb R^{m\times n}$05 and BUT and $A\in\mathbb R^{m\times n}$06 iterations without preconditioning (Li et al., 23 May 2025).

6. Backward error, conditioning, and constrained extensions

The backward-error theory of ILS centers on the normwise backward error

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

defined for a computed vector $A\in\mathbb R^{m\times n}$08. A tight linearization estimate is obtained from the Jacobian

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

namely

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

For sufficiently small true backward error, $A\in\mathbb R^{m\times n}$11. The same paper derives explicit mixed and componentwise condition numbers for $A\in\mathbb R^{m\times n}$12 via dual-norm techniques, and shows that their tight upper bounds can be estimated cheaply by the classical $A\in\mathbb R^{m\times n}$13-norm power method while using the QR–Cholesky solver (Diao et al., 2016).

Condition analysis has also been developed in a partial and structured form. For a linear functional $A\in\mathbb R^{m\times n}$14, the partial unified condition number is expressed through the Fréchet derivative matrix $A\in\mathbb R^{m\times n}$15, and by choosing specific norms and weights one recovers partial normwise, mixed, and componentwise condition numbers. When $A\in\mathbb R^{m\times n}$16 and $A\in\mathbb R^{m\times n}$17 belong to prescribed linear subspaces, the structured condition numbers satisfy inequalities such as $A\in\mathbb R^{m\times n}$18. Three estimators are proposed: a probabilistic spectral-norm estimator based on Lanczos bidiagonalization, a small-sample $A\in\mathbb R^{m\times n}$19-norm estimator, and a small-sample $A\in\mathbb R^{m\times n}$20-norm estimator. In the Toeplitz-structured experiments, the structured partial $A\in\mathbb R^{m\times n}$21-norm condition number could be $A\in\mathbb R^{m\times n}$22–$A\in\mathbb R^{m\times n}$23 smaller than the unstructured one (Li et al., 2016).

The equality-constrained indefinite least squares problem extends ILS to

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

or equivalently $A\in\mathbb R^{m\times n}$25 in alternative notation, under rank and definiteness assumptions on the constraint nullspace. Its solution is characterized by an augmented saddle-point system, and its sensitivity is described by projected condition numbers that unify normwise, mixed, and componentwise measures. For the $A\in\mathbb R^{m\times n}$26-norm case, the theory provides compact Kronecker-free formulas, while for mixed and componentwise quantities it provides tight upper bounds computable without Kronecker products (Wang et al., 2016). The normwise backward error for the constrained problem admits a linearized estimate $A\in\mathbb R^{m\times n}$27; if

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

then the exact backward error satisfies $A\in\mathbb R^{m\times n}$29, so the linearization is provably tight up to a small constant factor (Diao et al., 2018).

Across these developments, ILS is best understood not as a minor variant of least squares but as a distinct indefinite-signature problem class. Its normal equations resemble those of ordinary least squares, yet its algorithmic behavior is organized by signature-induced block structure, its fastest Krylov solvers depend on spectrum-clustering preconditioners, and its perturbation theory requires more care because naive normwise bounds can be unattainable (Grcar, 2010, Shaldehi et al., 28 Feb 2026).