Documents Live, a web authoring and publishing system
Documents Live, a web authoring and publishing system
If you see this, something is wrong
First published on Monday, Aug 24, 2026 and last modified on Monday, Aug 24, 2026
New York University Shanghai 567 West Yangsi Road Shanghai, 200126, P.R. China
We give a self-contained finite-precision analysis of the symmetric Lanczos algorithm without reorthogonalization. In particular, we derive the perturbed three-term recurrence, Paige’s loss-of-orthogonality identity, containment of all computed Ritz values, and localization of stabilized Ritz values. We then prove a Greenbaum-type backward stability result, exhibiting a nearby problem on which exact Lanczos produces the computed tridiagonal matrix. Our proofs simplify those of Paige and Greenbaum, at the cost of hiding polynomial factors in the iteration count.
Let \(A=A^\mathsf{T}\in\mathbb R^{n\times n}\) and let \( v\in\mathbb R^n\) be nonzero. In exact arithmetic, \(k\) non-terminating steps of the Lanczos algorithm (Algorithm 1) applied to \((A,v)\) produce a matrix \(V_k = [v_1, …, v_k]\) with orthonormal columns, \(v_1=v/\lVert v\rVert \), and a symmetric tridiagonal matrix \(T_k\in\mathbb R^{k\times k}\) satisfying
Hence, for \( j\leq k+1\) , the vectors \( {v_1, …, v_{j}}\) form an orthonormal basis for the Krylov subspace
The matrices \( V_k\) and \( T_k\) can then be used for downstream tasks such as estimating eigenvalues, solving linear systems of equations, and problems involving matrix functions [1, 2, 3, 4].
It was known since the introduction of the Lanczos algorithm that the behavior of the algorithm in finite precision arithmetic can be drastically different [5]. Indeed, the columns of \(V_k\) need not remain mutually orthogonal (and in fact can eventually become nearly linearly dependent), and the computed tridiagonal matrix may differ greatly from what would have been produced in exact arithmetic. Figure 1 illustrates this on a small example; the same run is revisited in Section 5, where the backward model is built for it explicitly.
Remarkably, this instability does not render the algorithm useless. Lanczos without reorthogonalization is used widely and successfully: Ritz values continue to converge to eigenvalues of \(A\) (although possibly with spurious copies), and approximations built from the computed \(T_k\), such as those to quadratic forms and matrix functions, are often nearly as accurate as their exact-arithmetic counterparts. Explaining this apparent paradox has occupied numerical analysts for decades. Much of our current understanding of the algorithm’s behavior in finite precision arithmetic is due to Chris Paige and Anne Greenbaum.
Paige showed that the algorithm is well-behaved locally, even in the presence of rounding errors [6, 7, 8]. In particular, he proved that the computed Lanczos vectors nearly satisfy a symmetric three-term recurrence, are nearly unit length, and that consecutive vectors are nearly orthogonal. He also proved that loss of orthogonality of the Lanczos vectors goes hand-in-hand with the convergence of a Ritz value (eigenvalue of \( T_k\) ) to an eigenvalue of \( A\) [9]; this result is often referred to as “Paige’s theorem”. Finally, it is commonly attributed to Paige that every computed Ritz value lies nearly within the interval containing the eigenvalues of \( A\) .
Greenbaum later proved a backward-stability-type result for Lanczos [10]. Building on the work of Paige, she showed that the computed tridiagonal matrix can also be obtained as the output of Lanczos run in exact arithmetic on a nearby problem instance. The nearby problem instance is of a larger dimension than \( n\) , but is close in the sense that its eigenvalues lie within small intervals about the eigenvalues of \( A\) , and the spectral weights of its starting vector on each such interval are close to those of the original starting vector [10, \S8]; see also [11].
Several resources exist for learning about these results. The original papers are rigorous and complete, but they are not an easy read: the analysis is spread across multiple works, the bounds are assembled through careful elementwise rounding-error bookkeeping, and Greenbaum herself describes the theorems and proofs of [10] as “gory”. The survey of Meurant and Strakoš [12] is an accessible guide to the results, their history, and their consequences, but it states the central theorems without proof, deferring to the original papers. Meurant’s book [1] collects these results and describes many aspects of their proofs carefully, but it likewise stops short of complete proofs of some key theorems. Parlett [13] distills a proof of Paige’s Theorem to its essential mechanism, but makes use of heuristic arguments to do so. In fact, in none of these resources do we find a proof of the Ritz value containment.
The aim of this note is to present a simplified, rigorous, self-contained analysis that captures the essence of the results of Paige and Greenbaum. Our main tools in this quest are two design choices:
While the resulting bounds are weaker than existing ones, we believe the tradeoff with simplicity is a good one. First, re-deriving standard floating-point bounds for linear-algebra primitives from the standard model of floating point arithmetic sheds little light on the behavior of Lanczos itself. Second, existing bounds for Lanczos are already numerically vacuous for even relatively small iteration counts \( k\) and dimensions \( n\) in double precision arithmetic, so their value is in the intuition they provide rather than as a usable bound. Third, polynomial dependencies on \( k\) do not significantly impact the theoretical number of bits of precision needed, which grows only logarithmically with the target accuracy.
Definition 1 (Polynomial convention)
Each occurrence of
is shorthand for \(x\le c k^dy\), where \(c,d>0\) are fixed absolute constants attached to that instance. The constant and exponent may change from one occurrence to the next.
Remark 1
This notation does not hide things that shouldn’t be hidden such as the dimension, spectral quantities of \( A\) , etc. In the companion Lean proof, available at https://github.com/tchen-research/simple-lanczos-stability , we provide explicit constants and exponents.
The remainder of this note is organized as follows. Section 2 states the arithmetic model and the computed recurrence. Section 3 derives our versions of Paige’s results: local error bounds, the loss-of-orthogonality theorem, and containment and localization of the computed Ritz values. Section 4 constructs the Greenbaum-type backward model.
We assume linear-algebra kernels can be performed to relative accuracy quantified by the accuracy parameter \( \varepsilon\) . The precise form of our assumptions is based on the bounds one obtains from the standard model of floating point arithmetic [14]. As such, these assumptions can be ensured by setting the machine precision sufficiently small relative to the ambient dimension \( n\) .
Assumption 1 (Arithmetic kernels)
For compatible vectors and scalars, the computed kernels satisfy
Copying and transposition are exact.
Assumption 2 (Normalization)
For \(y\ne\vec0\), the call \((\beta,v)\)=normalize\((y)\) returns \(\beta>0\) and \(v\) such that
If the vector to be normalized is zero, the algorithm terminates before making the call.
Finally, we assume the accuracy parameter \( \varepsilon\) is small.
Assumption 3 (Small \(\varepsilon\))
For a run of at most \( k_{\max}\) iterations, the accuracy parameter satisfies \( 0\leq \varepsilon \leq 1/(20k_{\max})\) .
We analyze Paige’s computational variant A1 [8], stated explicitly in Algorithm 1. This is the now-standard implementation of Lanczos.
Remark 2
The initial call returns \(v_1\) proportional to \(v+\vec\Delta_0\) with \(\lVert \vec\Delta_0\rVert \le\varepsilon\lVert v\rVert \), by the first equation in Assumption 2. In particular \(v_1\) satisfies the second equation in Assumption 2, like every later Lanczos vector, and this is the only property of \(v_1\) the analysis of Section 3 uses.
Fix a run of \(k\) completed iterations for which \(\beta_1,…,\beta_{k-1}>0\); the last step may be nonterminal (\(\beta_k>0\)) or terminal (\(\beta_k=0\)). Let
Thus \(T_k=T_k^\mathsf{T}\), and its internal subdiagonals are positive. If the run terminates at step \(k\), then \(z_k=\vec0\) and \(\beta_k=0\). In that case we adopt the conventions \(\vec\Delta_k=\vec0\) and \(\beta_kv_{k+1}=\vec0\), even though \(v_{k+1}\) is not formed.
Define the exact errors of the five kernel calls and record their bounds in one panel; the normalization call carries two guarantees, so it contributes the final two lines:
(1)
(2)
(3)
(4)
(5)
(6)
These are identities together with their consequences from the arithmetic model. At a terminal final step, Equation 5 holds by the convention above rather than by a normalization call, and Equation 6 does not apply since \(v_{k+1}\) is not formed.
We begin by showing that the Lanczos algorithm behaves well locally; i.e. the computed Lanczos vectors nearly satisfy a symmetric three-term recurrence, are nearly unit length, and that consecutive vectors are nearly orthogonal. Such bounds appeared in Paige’s early work; see e.g. [8, \S2-3]. Specifically, we bound the following quantities:
(7)
Theorem 1 (Local error bounds)
For \(1\le j\le k\),
The bound on \(g_j\) is immediate from Assumption 2. The other two are not: \(f_j\) aggregates the kernel errors of step \(j\), while \(p_j\) couples the vectors \(v_j\) and \(v_{j+1}\), which come from different normalization calls, scaled by \(\beta_j\), which can be as large as \(\lVert A\rVert \). Both bounds rest on an a priori envelope for the computed quantities; we establish the envelope first and prove Theorem 1 after it.
Lemma 1 (A priori coefficient bounds)
For \(1\le j\le k\),
Proof
The two computed projection equations in the local-error panel give
The bracketed rank-one update does not increase the norm, even though \(v_j\) is only nearly normalized. By Assumption 2, \(\lVert v_j\rVert ^2\le1+\varepsilon<2\). Therefore, direct expansion gives
The inner-product model and Cauchy–Schwarz give
The update error therefore satisfies
Combining these three bounds gives
where the last inequality holds because \(\varepsilon\le1/20\), by Assumption 3.
It remains to control the size of \(w_j\). The first two rows of the local-error panel and \(\lVert v_i\rVert \le\sqrt{1+\varepsilon}\) give
and
At a nonterminal step, Equation 5 and \(\lVert v_{j+1}\rVert \ge\sqrt{1-\varepsilon}\), from Equation 6, give
again using \(\varepsilon\le1/20\). With the bound on \(z_j\), this yields the scalar recurrence
where the last line uses the two elementary inequalities
both valid for \(0\le\varepsilon\le1/20\). Starting from \(\beta_0=0\), induction on \(j\) then gives \(\beta_j\le4j\lVert A\rVert \): assuming the bound at \(j-1\),
because \(40(j-1)\varepsilon\le40k_{\max}\varepsilon\le2\) by Assumption 3. This is the only place the smallness assumption is used quantitatively: each step adds at most \(2\lVert A\rVert \) and inflates the previous bound by a factor \(1+10\varepsilon\), and the threshold \(\varepsilon\le1/(20k_{\max})\) is exactly what keeps the drift compounded over \(k_{\max}\) steps below the additive increment. Substituting \(\beta_{j-1}\le4(j-1)\lVert A\rVert \) into the displayed bounds, and again using \(\varepsilon\le1/20\), gives
Every quantity in the statement is therefore at most \(7j\lVert A\rVert \), which proves the claimed bounds. At a terminal final step \(\beta_k=0\), and the same bounds apply using the already controlled \(\beta_{k-1}\).
Proof (Proof of Theorem 1)
The bound on \(g_j\) is the unit-norm guarantee of the normalization calls: \(g_j=\delta_j\) in Equation 6 for \(j\ge2\), and the initial call gives the same bound for \(g_1\).
For \(f_j\), eliminate the intermediate vectors: substitute Equation 1 into Equation 2, the result into Equation 4, and that into Equation 5. Rearranging the result against the definition of \(f_j\) in Equation 7 gives
Apply the triangle inequality and insert the four error bounds from the local-error panel:
By Assumption 2 each \(\lVert v_i\rVert \) is at most \(\sqrt{1+\varepsilon}\), and by Lemma 1 each of \(\lVert u_j\rVert \), \(\lVert w_j\rVert \), \(\lvert \alpha_j\rvert \), \(\lVert z_j\rVert \), and \(\beta_{j-1}\) is \(\lesssim_k\lVert A\rVert \). The bracket is therefore \(\lesssim_k\lVert A\rVert \), which proves the bound on \(f_j\).
For \(p_j\), \crefrange{eq:local-alpha}{eq:local-normalize} give
Since \(v_j^\mathsf{T}w_j=\alpha_j-e^\alpha_j\), the two exact copies of \(\alpha_j\) cancel, leaving
Apply Assumption 2, Equation 3, Equation 4, Equation 5, and Lemma 1 to obtain the bound on \(p_j\). The matrix–vector error has already disappeared in the exact cancellation.
Stack the recurrence errors into \(F_k:=[f_1,…,f_k]\) and expresses the run as a perturbed block recurrence,
(8)
Define the Gram defect and its triangular decomposition by
(9)
where \(R_k:=\operatorname{up}(\vec\Omega_k)\) and \(D_k:=\operatorname{diag}(g_1,…,g_k)\), and \(\operatorname{up}\) retains the strict upper triangle. No smallness is assumed for \(R_k\) or \(\vec\Omega_k\). Set
(10)
The following identity, which prevents the Gram defect \(\vec\Omega_k\) from being arbitrary by forcing it to nearly commute with \(T_k\), is the algebraic core of Paige’s analysis [8, eq. (22)], [9, eqs. (2.17),(2.18)].
Lemma 2 (Gram commutator)
The governing relation implies the exact identity
Every entry of \(G_k\) has magnitude \(\lesssim_k\varepsilon\lVert A\rVert \). Moreover, with \(p_0=0\) and
the corresponding upper-triangular identity is
The last term is supported only on the first superdiagonal, and every entry of it also has magnitude \(\lesssim_k\varepsilon\lVert A\rVert \).
Proof
Premultiplying Equation 8 by \(V_k^\mathsf{T}\) gives
The left side is symmetric. Equating the right side with its transpose proves the first identity. The entry bound follows from \((G_k)_{rs}=v_r^\mathsf{T}f_s-f_r^\mathsf{T}v_s\), Cauchy–Schwarz, and Theorem 1, which controls both factors: \(\lVert v_r\rVert \le\sqrt{1+\varepsilon}\) through the bound on \(g_r\), and \(\lVert f_s\rVert \lesssim_k\varepsilon\lVert A\rVert \).
Insert Equation 9 in the first identity and retain the strict upper triangle. Tridiagonality of \(T_k\) and strict upper triangularity of \(R_k\) imply that \(T_kR_k^\mathsf{T}-R_k^\mathsf{T} T_k\) is lower triangular. The same two properties pin each diagonal entry of \(T_kR_k\) and \(R_kT_k\) to a single product:
where the first is zero for \(j=1\) and the second is absent for \(j=k\), because \(R_k\) has no column \(k+1\). Thus the diagonal of \(T_kR_k-R_kT_k\) is \(p_{j-1}-p_j\), except that its last entry is \(p_{k-1}\); the right side of the second identity agrees, because \(N_k\) contributes \(p_{j-1}-p_j\) and \((C_k)_{kk}=p_k\) restores the missing \(-p_k\). This gives the second identity. Finally, \((T_kD_k-D_kT_k)_{j,j+1}=\beta_j(g_{j+1}-g_j)\), whose magnitude is at most \(2\varepsilon\beta_j\lesssim_k\varepsilon\lVert A\rVert \), by Theorem 1, Lemma 1.
Let \((\theta,y)\) be a unit eigenpair of \(T_k\), and define the corresponding Ritz vector \(x:=V_ky\). Because \(T_ky=\theta y\), the governing relation Equation 8 gives the exact residual formula
(11)
and therefore
(12)
In exact arithmetic \(F_k=\vec0\) and \(\lVert v_{k+1}\rVert =1\), so \(\beta_k\lvert e_k^\mathsf{T}y\rvert \) is exactly the residual norm of the pair \((\theta,x)\); in finite precision it determines the residual norm up to \(\lesssim_k\varepsilon\lVert A\rVert \), while remaining computable from \(T_k\) alone. We therefore call \(\beta_k\lvert e_k^\mathsf{T}y\rvert \) the scalar residual estimate of the pair. When the residual estimate is small, the Ritz pair is nearly an eigenpair: since \(A\) is symmetric, Equation 12 gives
so \(\theta\) is then nearly an eigenvalue of \(A\), and \(x/\lVert x\rVert \) is nearly a corresponding eigenvector whenever that eigenvalue is well separated from the rest of the spectrum. The caveat is the denominator: loss of orthogonality can make \(\lVert x\rVert \) small, and removing this caveat is the purpose of the next subsection. The next theorem is the central loss-of-orthogonality bound in Paige’s finite-precision analysis [9, eqs. (3.11)–(3.13)].
Theorem 2 (Paige loss-of-orthogonality bound)
For every unit eigenpair \((\theta,y)\) of \(T_k\), the Ritz vector \(x=V_ky\) satisfies
At a nonterminal step, Theorem 2 gives the usual convergence–orthogonality dichotomy:
(13)
Thus appreciable loss of orthogonality against a Ritz direction can occur only when its scalar residual estimate is small, that is, by Equation 12, only when the Ritz pair has nearly converged. The same tradeoff is visible numerically in Figure 2.
Proof
Fix a unit eigenpair \((\theta,y)\) and sandwich the second identity in Lemma 2 between \(y^\mathsf{T}\) and \(y\). The left side vanishes because \(T_ky=\theta y\). The \(C_k\) term is the boundary product,
so the sandwiched identity rearranges to the exact expression
Because \(y\) is a unit vector, each quadratic form on the right is bounded by the operator norm, and hence by the Frobenius norm, of its matrix. All three matrices are \(k\times k\) with entries of magnitude \(\lesssim_k\varepsilon\lVert A\rVert \): the entries of \(N_k\) are differences of the local products \(p_j\), which Theorem 1 controls, and the other two are covered by Lemma 2. Each Frobenius norm is therefore \(\lesssim_k\varepsilon\lVert A\rVert \), and the theorem follows.
For every \(1\le t\le k\), the tridiagonal \(T_t\) produced after the first \(t\) iterations satisfies the preceding identities and bounds; in particular, Theorem 2 applies to the unit eigenpairs of every \(T_t\). Paige analyzes some of the consequences of this observation in [9, \S3].
Lemma 3 (Paige descent)
Let \((\theta,y)\) be a unit eigenpair of \(T_t\). If \(\lvert y^\mathsf{T} R_ty\rvert >3/8\), then there are \(1\le r\le s<t\) and a unit eigenpair \((\theta_r^{(s)},y_r^{(s)})\) of \(T_s\) such that
Proof
Put \(\rho:=y^\mathsf{T} R_ty\). For \(s<t\), write \(y^{[s]}\) for the first \(s\) entries of \(y\), and expand it in an orthonormal eigenbasis \(y_1^{(s)},…,y_s^{(s)}\) of \(T_s\): \(y^{[s]}=\sum_{r=1}^{s} c_r^{(s)}y_r^{(s)}\). The leading principal \(s\times s\) block of \(T_t\) is \(T_s\), and, by tridiagonality, the only entry of the first \(s\) rows of \(T_t\) outside this block is \(\beta_s\), in position \((s,s+1)\). The first \(s\) rows of \(T_ty=\theta y\) therefore read
Multiplying on the left by \((y_r^{(s)})^\mathsf{T}\), and using \((y_r^{(s)})^\mathsf{T}T_s=\theta_r^{(s)}(y_r^{(s)})^\mathsf{T}\) together with \((y_r^{(s)})^\mathsf{T}y^{[s]}=c_r^{(s)}\), gives the division-free identity
(14)
Since \(R_t\) is the strict upper triangle of \(V_t^\mathsf{T}V_t-I\), its entries are \((R_t)_{ij}=v_i^\mathsf{T}v_j\) for \(i<j\). Grouping the quadratic form \(\rho\) by the column index \(j\) gives
where the last equality reindexes by \(s=j-1\) and uses \(\sum_{i\le s}(e_i^\mathsf{T}y)v_i=V_sy^{[s]}\). Substituting the eigenbasis expansion of \(y^{[s]}\) splits each term over the level-\(s\) Ritz vectors \(V_sy_r^{(s)}\):
(15)
This expansion is [9, eqs. (3.18),(3.19)]. Each factor \((V_sy_r^{(s)})^\mathsf{T}v_{s+1}\) is the orthogonality defect of the Lanczos vector \(v_{s+1}\) against a Ritz vector of an earlier iteration—precisely the quantity that Theorem 2 couples to the corresponding residual estimate. Introduce the Paige products and their maximum
(16)
Each \(P_{rs}\) is the boundary product of the unit eigenpair \((\theta_r^{(s)},y_r^{(s)})\) of \(T_s\), so Theorem 2 applied to the first \(s\) iterations gives \(\mathcal E_*\lesssim_k\varepsilon\lVert A\rVert \). The products \(\tau_{rs}\) and \(P_{rs}\) share the factor \((V_sy_r^{(s)})^\mathsf{T}v_{s+1}\): multiplying \(\tau_{rs}\) by \(\theta-\theta_r^{(s)}\) and applying Equation 14, or multiplying by \(\beta_s(e_s^\mathsf{T}y_r^{(s)})\) and rearranging factors, gives the exact relations
(17)
Set \(\Psi:=8t\,\mathcal E_*\lesssim_k\varepsilon\lVert A\rVert \). We prove the conclusion of the lemma, with both right sides replaced by \(\Psi\), in contrapositive form: assuming that every pair \((r,s)\) satisfies
we show \(\lvert \rho\rvert <3/8\), so that the hypothesis of the lemma fails.
If \(\mathcal E_*=0\), then \(\Psi=0\) and every \(P_{rs}\) vanishes; for each pair, dividing the relation in Equation 17 that corresponds to the assumed strict inequality by its nonzero left factor gives \(\tau_{rs}=0\), so \(\rho=0\) by Equation 15. We may therefore assume \(\mathcal E_*>0\). For each summand in Equation 15, either \(\lvert \theta-\theta_r^{(s)}\rvert >\Psi\), in which case the first relation in Equation 17 gives
or \(\beta_s\lvert e_s^\mathsf{T}y_r^{(s)}\rvert >\Psi\), in which case the second relation in Equation 17 gives
Therefore
The two sums are bounded separately, freely using \(s\le t\) and \(\lVert y^{[s]}\rVert \le\lVert y\rVert =1\). First,
Second, by Cauchy–Schwarz over \(r\) and orthonormality of the eigenbasis,
so Cauchy–Schwarz over \(s\) gives
Since \(\mathcal E_*/\Psi=1/(8t)\),
In both cases \(\lvert \rho\rvert <3/8\), which completes the contrapositive.
Remark 3
Paige reaches the same conclusion by introducing the ratio \((e_{s+1}^\mathsf{T}y)/(\theta-\theta_r^{(s)})\) and bounding a Frobenius norm [9, eqs. (3.30)–(3.34)]. Splitting instead on which of \(\lvert \theta-\theta_r^{(s)}\rvert \) and \(\beta_s\lvert e_s^\mathsf{T}y_r^{(s)}\rvert \) is large keeps every step division-free.
The next theorem states that every computed Ritz value lies in a slightly enlarged spectral interval of \(A\). This fact is essential to the analysis of Lanczos-based methods for matrix functions [15, 16, 17]; see [4].
Theorem 3 (Containment of every computed Ritz value)
Every \(\theta\in\operatorname{spec}(T_k)\) satisfies
Remark 4
Paige writes this bound [9, Eq. 3.48], but does not include a proof. Many papers and books reference [9, Eq. 3.48] but also do not provide the proof. In fact, the only proof of this result of which we are aware of is [18, Theorem A.1] via an entirely different argument.
Proof
Let \(T_ky=\theta y\), \(\lVert y\rVert =1\), put \(x=V_ky\), and set \(\rho=y^\mathsf{T} R_ky\). Suppose first that \(\lvert \rho\rvert \le3/8\); we show the Ritz vector cannot then be too short. Expanding \(V_k^\mathsf{T}V_k\) via the Gram splitting Equation 9, and using \(y^\mathsf{T}R_k^\mathsf{T}y=y^\mathsf{T}R_ky=\rho\),
The diagonal term satisfies \(\lvert y^\mathsf{T}D_ky\rvert \le\max_{j\le k}\lvert g_j\rvert \le\varepsilon\) by Theorem 1, so
using \(\varepsilon\le1/20\), as ensured by Assumption 3. Premultiply Equation 11 by \(x^\mathsf{T}\). The Rayleigh quotient \(x^\mathsf{T} Ax/\lVert x\rVert ^2\) lies in the spectral interval of \(A\), so
by Theorem 2 and \(\lVert F_k\rVert \le\sqrt{k}\,\max_{j\le k}\lVert f_j\rVert \lesssim_k\varepsilon\lVert A\rVert \) from Theorem 1.
If \(\lvert \rho\rvert >3/8\), we descend to an earlier iteration at which the first branch applies. Set \(t_0:=k\) and \((\theta^{(t_0)},y^{(t_0)}):=(\theta,y)\), and iterate the following step. Given a unit eigenpair \((\theta^{(t_i)},y^{(t_i)})\) of \(T_{t_i}\) with \(\lvert (y^{(t_i)})^\mathsf{T}R_{t_i}y^{(t_i)}\rvert >3/8\), Lemma 3 applied to \(T_{t_i}\) produces an index \(t_{i+1}<t_i\) and a unit eigenpair \((\theta^{(t_{i+1})},y^{(t_{i+1})})\) of \(T_{t_{i+1}}\) with
The process stops at the first index \(m\) with \(\lvert (y^{(t_m)})^\mathsf{T}R_{t_m}y^{(t_m)}\rvert \le3/8\). It must stop: the indices \(k=t_0>t_1>…\ge1\) strictly decrease, and \(R_1\), the strict upper triangle of a \(1\times1\) matrix, is zero, so the stopping condition holds at \(t=1\) at the latest; in particular, \(m\le k-1\).
Every displacement above is an instance of the same bound in Lemma 3, so all carry the same absolute constants \(c\) and \(d\): each step satisfies \(\lvert \theta^{(t_i)}-\theta^{(t_{i+1})}\rvert \le ck^d\varepsilon\lVert A\rVert \). The triangle inequality over the at most \(k-1\) steps therefore gives
The stopping pair satisfies the hypothesis of the first branch, applied to the first \(t_m\) iterations, so \(\operatorname{dist}(\theta^{(t_m)}, [\lambda_{\min}(A),\lambda_{\max}(A)]) \lesssim_k\varepsilon\lVert A\rVert \). Adding the two displacements proves the theorem.
Theorem 3 places every computed Ritz value near the spectral interval of \(A\). The final theorem of this section sharpens this to the spectrum itself for Ritz values whose scalar residual estimate \(\beta_k\lvert e_k^\mathsf{T}y\rvert \) is small (in Paige’s terminology, Ritz values that have stabilized) [9]. In exact arithmetic such a bound would follow at once from Equation 12; in finite precision that route divides by \(\lVert x\rVert \), which loss of orthogonality can make small. Paige’s treatment of stabilized and clustered Ritz values was subsequently refined by Wülling [19]; we do not pursue such refinements here.
Theorem 4 (Localization of stabilized Ritz values)
Let \((\theta,y)\) be a unit eigenpair of \(T_k\). Then
In particular, if the recurrence terminates at step \(k\), so \(\beta_k=0\), every eigenvalue of \(T_k\) satisfies \(\operatorname{dist}(\theta,\operatorname{spec}(A))\lesssim_k\varepsilon\lVert A\rVert \).
Proof
Put \(r:=\beta_k\lvert e_k^\mathsf{T}y\rvert \). Starting from \((\theta,y)\), apply Lemma 3 whenever \(\lvert y^\mathsf{T} R_ty\rvert >3/8\). The iteration number strictly decreases, so after at most \(k-1\) steps the process reaches an eigenpair \((\widehat\theta,\widehat{y})\) of some \(T_t\) whose Ritz vector \(\widehat{x}=V_t\widehat{y}\) has \(\lVert \widehat{x}\rVert ^2\ge1/5\). As in the proof of Theorem 3, the accumulated displacement satisfies \(\lvert \theta-\widehat\theta\rvert \lesssim_k\varepsilon\lVert A\rVert \). The residual estimate of the final pair is \(r\) if no descent occurred, and otherwise is \(\lesssim_k\varepsilon\lVert A\rVert \) by Lemma 3; in either case it is \(\lesssim_k r+\varepsilon\lVert A\rVert \). Since \(F_t\) consists of the first \(t\) columns of \(F_k\), \(\lVert F_t\rVert \le\lVert F_k\rVert \le\sqrt{k}\,\max_{j\le k}\lVert f_j\rVert \lesssim_k\varepsilon\lVert A\rVert \). Applying Equation 12 to the final pair gives
Adding the accumulated displacement proves the theorem; at termination \(r=0\).
We now show that the computed \(T_k\) is the exact output of \(k\) steps of Lanczos applied to a symmetric matrix that is close to an orthogonal copy of \(I_k\otimes A\), a backward interpretation of the kind introduced by Greenbaum in [10]. Our construction differs from Greenbaum’s. Whereas Greenbaum builds the model by extending the computed tridiagonal matrix, continuing the recurrence in exact arithmetic and tracking where the extended eigenvalues fall [10, §\S5–7], we instead orthogonalize the computed Lanczos vectors exactly, by a dilation into \(k\) stacked copies of \(\mathbb R^n\), and then complete the dilated block to a symmetric matrix. The dilation has three payoffs: the model is built in closed form, in dimension exactly \(nk\), with no recurrence to continue; its perturbation is of the same order \(\varepsilon\lVert A\rVert \) as the local error bounds (Greenbaum’s bound depends on the fourth-root of the accuracy parameter [10, eq. (7.10)]); and the starting spectral weights are preserved exactly, so that all error is concentrated in a single perturbation bound. Moreover, the construction does not assume that \(V_k\) is well-conditioned, that \(k\) precedes exact breakdown for \((A,v_1)\), or that the spectral measure is regular.
The construction takes as input the governing recurrence in Equation 8 together with the local error bounds of Theorem 1; nothing else about the run is used. It proceeds in three steps that mirror the structure of Section 3. First, we rescale the computed Lanczos vectors to have exactly unit norm and check that the local error bounds survive, in parallel with Theorem 1. Second, we bound a commutator built from the strict upper triangle of the normalized Gram matrix, in parallel with Lemma 2. Third, we orthogonalize exactly: the normalized vectors and the Gram triangle assemble into an isometry into a space of \(k\) stacked copies of \(\mathbb R^n\), on which the recurrence is reproduced up to a residual of size \(\lesssim_k\varepsilon\lVert A\rVert \).
Throughout this subsection, suppose Algorithm 1 completes \(k\) iterations under the assumptions of Section 2. Let \(d_j:=\lVert v_j\rVert \), \(D:=\operatorname{diag}(d_1,…,d_k)\), and \(\bar{V}:=V_kD^{-1}\), so the columns \(\bar{v}_j=v_j/d_j\) of \(\bar{V}\) are exactly unit vectors. Write \(\beta_kv_{k+1}/d_k=\bar\beta_k\bar{v}_{k+1}\), where \(\bar\beta_k\ge0\) and \(\bar{v}_{k+1}\) is a unit vector, chosen arbitrarily if \(\beta_kv_{k+1}=\vec0\). In this basis the governing recurrence in Equation 8 reads
(18)
which defines the normalized residual \(E\). Finally, in parallel with Equation 9, decompose the normalized Gram matrix as
(19)
where \(U\) is strictly upper triangular; the diagonal is exactly the identity because the columns are exactly unit.
The first step is the analogue of Theorem 1 for the normalized quantities.
Lemma 4 (Normalized local errors)
The normalized residual satisfies \(\lVert E\rVert \lesssim_k\varepsilon\lVert A\rVert \), and the normalized weighted adjacent products satisfy
Proof
Since \(d_j^2=1+g_j\), the bound \(\lvert g_j\rvert \le\varepsilon\) of Theorem 1 gives \(\lVert D-I\rVert \le\varepsilon\); moreover, \(\varepsilon\le3/4\) by Assumption 3, so \(d_j^2\ge1-\varepsilon\ge1/4\) and \(\lVert D^{-1}\rVert \le2\). Combining the two, \(\lVert D^{-1}-I\rVert =\lVert D^{-1}(I-D)\rVert \le2\varepsilon\). Substituting \(V_k=\bar{V}D\) into Equation 8 and right multiplying by \(D^{-1}\) identifies the residual in Equation 18 as
Since \(\lVert T_k\rVert \lesssim_k\lVert A\rVert \) by Lemma 1,
while \(\lVert F_kD^{-1}\rVert \le2\sqrt{k}\,\max_{1\le j\le k}\lVert f_j\rVert \lesssim_k\varepsilon\lVert A\rVert \) by Theorem 1. Since \(\bar{V}\) has unit columns, \(\lVert \bar{V}\rVert \le\sqrt{k}\), and hence, combined with the two preceding bounds, \(\lVert E\rVert \lesssim_k\varepsilon\lVert A\rVert \).
Up to the normalization scalars, each adjacent product is a weighted adjacent product from Equation 7,
and \(d_j\ge1/2\), so both are \(\lesssim_k\varepsilon\lVert A\rVert \) by Theorem 1.
The second step bounds a commutator built from \(U\), in parallel with the upper-triangular identity in Lemma 2: there the strict upper triangle \(R_k\) of the Gram defect satisfied a commutator relation driven by the local errors, and here \(U\) plays the role of \(R_k\).
Lemma 5 (Normalized Gram commutator)
The matrix \(T_kU-UT_k-\bar\beta_kae_k^\mathsf{T}\) is upper triangular, and its norm is \(\lesssim_k\varepsilon\lVert A\rVert \).
Proof
Put \(M:=T_kU-UT_k-\bar\beta_kae_k^\mathsf{T}\). Tridiagonality of \(T_k\) and strict upper triangularity of \(U\) make \(M\) upper triangular. Left multiplication of Equation 18 by \(\bar{V}^\mathsf{T}\) gives
The left side is symmetric, so the right side equals its own transpose:
Substituting \(G=I+U+U^\mathsf{T}\), cancelling the \(T_k\) terms coming from the identity, and collecting the remaining terms into \(M\) and \(M^\mathsf{T}\) gives
(20)
which controls the strict upper triangle of \(M\). Its diagonal entries are
with absent endpoint terms omitted; each is \(\lesssim_k\varepsilon\lVert A\rVert \) by the adjacent-product bounds of Lemma 4. Since \(M^\mathsf{T}\) is lower triangular, every strict upper entry of \(M\) equals the corresponding entry of \(M-M^\mathsf{T}\), and reading that entry off Equation 20 gives, for \(i<j\),
because \(\bar{v}_i\) and \(e_j\) are unit vectors, and \(\lVert E\rVert \lesssim_k\varepsilon\lVert A\rVert \) by Lemma 4. Since \(M\) is \(k\times k\) and upper triangular, these entrywise bounds give \(\lVert M\rVert \lesssim_k\varepsilon\lVert A\rVert \).
The third step contains no error analysis at all; it is exact linear algebra, using only that \(\bar{V}\) has unit columns and \(\bar{v}_{k+1}\) is a unit vector. Where Paige used the triangle \(R_k\) of the Gram defect to track the loss of orthogonality, we use the triangle \(U\) to undo it: setting
(21)
the recombination \(\bar{V}C\) of the columns of \(\bar{V}\) fails to be an isometry only by the auxiliary block \(K\).
Lemma 6 (Nilpotent isometry)
The matrix \(K\) is strictly upper triangular, so \(K^k=\vec0\); moreover \(\lVert K\rVert \le1\) and \(\lVert C\rVert \le2\). In addition,
Proof
Since \(K=UC=CU\), it is strictly upper triangular, and hence \(K^k=\vec0\). Since \(U=C^{-1}-I\), the Gram decomposition in Equation 19 reads
Multiplying by \(C^\mathsf{T}\) on the left and \(C\) on the right, the inverses cancel term by term:
Expanding \(K^\mathsf{T}K=(I-C)^\mathsf{T}(I-C)=I-C-C^\mathsf{T}+C^\mathsf{T}C\) identifies the right side, so
(22)
Since \((\bar{V}C)^\mathsf{T}(\bar{V}C)=C^\mathsf{T}GC\), this is the first identity of the lemma. Since the Gram matrix \(G\) is positive semidefinite, \(K^\mathsf{T}K=I-C^\mathsf{T}GC\preceq I\), so \(\lVert K\rVert \le1\) and \(\lVert C\rVert =\lVert I-K\rVert \le2\); no such a priori bound holds for \(U\) itself.
For the second identity,
where the last equality uses Equation 22 and \(C+K=I\). Finally,
because \(C^\mathsf{T}(G+I)C=C+C^\mathsf{T}\).
With the three steps in place, the dilation is assembled by stacking powers of \(K\).
Lemma 7 (Physical dilation)
On \(k\) copies of \(\mathbb R^n\), let \(\widehat{A}=I_k\otimes A\), and write \(\operatorname{col}(…)\) for the map that stacks its arguments into a block column. There exist \(\widehat{V}\in\mathbb R^{nk\times k}\), \(\widehat{v}\in\mathbb R^{nk}\), and \(\widehat{E}\in\mathbb R^{nk\times k}\) such that, with \(\bar\beta_k\ge0\) as above,
Proof
Define
with \(\widehat{v}=y\) when \(k=1\). Since block columns multiply block by block, the identities of Lemma 6 turn every inner product involving \(\widehat{V}\) and \(\widehat{v}\) into a telescoping sum. First,
using \(K^k=\vec0\) in the last step. Similarly,
and
Also \(Ce_1=e_1\) and \(Ke_1=\vec0\), so \(\widehat{V}e_1=\operatorname{col}(\bar{v}_1,\vec0,…,\vec0)\), proving the second equation of the lemma.
It remains to propagate the recurrence. Write \(\widehat{V}_r:=\bar{V}CK^r\) and \(\widehat{v}_r\) for the \(r\)th blocks of \(\widehat{V}\) and \(\widehat{v}\), indexed from \(r=0\), so that \(\widehat{v}_0=y\) and \(\widehat{v}_r=\widehat{V}_{r-1}s\) for \(r\ge1\). Put \(M:=T_kU-UT_k-\bar\beta_kae_k^\mathsf{T}\) and
(23)
where the norm bound combines \(\lVert M\rVert \lesssim_k\varepsilon\lVert A\rVert \) from Lemma 5 with \(\lVert C\rVert \le2\) from Lemma 6. Rearranged, Equation 23 is an exchange rule: moving \(K\) leftward past \(T_k\) costs a boundary term and a small residual,
This is the only commuting we need. We claim that, for \(0\le r<k\),
(24)
stacking these \(k\) block rows gives the third equation of the lemma, with \(\lVert \widehat{E}\rVert \le\sqrt{k}\,\max_r\lVert E_r\rVert \lesssim_k\varepsilon\lVert A\rVert \).
The claim follows by induction on \(r\). For \(r=0\), right multiply Equation 18 by \(C\): since \(e_k^\mathsf{T}C=e_k^\mathsf{T}\), and \(T_kC=CT_k-\bar\beta_kse_k^\mathsf{T}-L\) by the exchange rule and \(C=I-K\),
and \(\bar{v}_{k+1}-\bar{V}s=y=\widehat{v}_0\), while \(\lVert E_0\rVert \le2\lVert E\rVert +\sqrt{k}\,\lVert L\rVert \lesssim_k\varepsilon\lVert A\rVert \) by Lemma 4. For the step from \(r\) to \(r+1\), right multiply Equation 24 by \(K\). The boundary term dies because \(e_k^\mathsf{T}K=\vec0\) (the last row of a strictly upper triangular matrix is zero), and the exchange rule moves the trailing \(K\) past \(T_k\), recreating the next boundary:
which is Equation 24 for \(r+1\), with \(\widehat{V}_rs=\widehat{v}_{r+1}\) and \(E_{r+1}=E_rK+\widehat{V}_rL\). Since \(\lVert K\rVert \le1\) and \(\lVert \widehat{V}_r\rVert \le\lVert \bar{V}\rVert \lVert C\rVert \lVert K\rVert ^r\le2\sqrt{k}\), each step grows the residual by at most \(2\sqrt{k}\,\lVert L\rVert \), so \(\lVert E_r\rVert \le\lVert E_0\rVert +2r\sqrt{k}\,\lVert L\rVert \lesssim_k\varepsilon\lVert A\rVert \).
Lemma 8 (Completion of an approximate initial Lanczos block)
Let \(\widehat{V}\), \(\widehat{v}\), and \(\widehat{E}\) be as in Lemma 7, and put \(N=nk\). There are an orthogonal matrix \(\widehat{U}\in\mathbb R^{N\times N}\) whose first \(k\) columns are the columns of \(\widehat{V}\), and a symmetric tridiagonal matrix \(\widetilde{T}\in\mathbb R^{N\times N}\), such that
The matrix \(\widetilde{T}\) may be reducible only after its leading \(k\times k\) block.
Proof
From the first and third equations of Lemma 7, \(\widehat{V}^\mathsf{T}\widehat{E}=\widehat{V}^\mathsf{T}\widehat{A}\widehat{V}-T_k\), which is symmetric. Hence
(25)
is symmetric, satisfies \(\vec\Delta\widehat{V}=-\widehat{E}\), and obeys \(\lVert \vec\Delta\rVert \le3\lVert \widehat{E}\rVert \). Thus \(B=\widehat{A}+\vec\Delta\) satisfies the exact recurrence
If \(\bar\beta_k\lVert \widehat{v}\rVert >0\), normalize \(\widehat{v}\) and absorb its norm into the next subdiagonal; if \(\bar\beta_k\lVert \widehat{v}\rVert =0\), \(\operatorname{range}(\widehat{V})\) is invariant. Continue exact Lanczos for \(B\). At any later breakdown, symmetry makes the orthogonal complement of the generated invariant subspace invariant, so tridiagonalize that complement independently. This gives an orthogonal \(\widehat{U}\), with first \(k\) columns \(\widehat{V}\), for which \(\widetilde{T}=\widehat{U}^\mathsf{T} B\widehat{U}\) is tridiagonal and has leading block \(T_k\). Finally,
Theorem 5 (Greenbaum-type backward model)
Suppose Algorithm 1 completes \(k\) iterations under the assumptions of Section 2. Put \(N=nk\). There exist symmetric \(N\times N\) matrices \(S\) and \(\widetilde{T}\) such that
Consequently, exact \(k\)-step Lanczos applied to \((\widetilde{T},e_1)\) returns the computed matrix \(T_k\) exactly.
Proof
Apply Lemma 7 and then Lemma 8. Let \(\widehat{U}\) be the resulting orthogonal matrix and set
The completion bound and the residual bound in Lemma 7 prove (iii). Moreover, \(\widehat{U}e_1=\widehat{V}e_1 =\operatorname{col}(\bar{v}_1,\vec0,…,\vec0)\) by the second equation in Lemma 7, proving (ii). Since the run completed \(k\) iterations, the internal subdiagonals \(\beta_1,…,\beta_{k-1}\) of \(T_k\) are positive, so exact Lanczos from \(e_1\) reads the leading block of \(\widetilde{T}\) through step \(k\).
This is a backward interpretation at the level of the starting spectral measure. It does not assert the componentwise spectral-weight comparison for every computed Lanczos vector given in [10, \S8]. Its advantages are that the perturbation is of the same order \(\varepsilon\lVert A\rVert \) as the local errors and that no inverse subdiagonal, starting weight, or spectral gap enters the construction. Our backwards matrix, which is nearly the direct sum of \( A\) with itself many times, is intuitively similar to the models used in [11] to “predict” the behavior of finite precision Lanczos.
We illustrate the backward model of Section 4 on the run shown in Figure 1. The matrix \(A\) is the \(n=64\) Strakoš matrix [20]: diagonal, with eigenvalues (in decreasing order) \(\lambda_i=\lambda_n+\frac{n-i}{n-1}(\lambda_1-\lambda_n)\rho^{i-1}\), where \(\lambda_1=1\), \(\lambda_n=10^{-3}\), and \(\rho=0.9\). Algorithm 1 is run for \(k=30\) iterations in IEEE single precision (\(\varepsilon\approx10^{-7}\)) from a random starting vector, and the backward model is then assembled in double precision, following the construction of Section 4 verbatim. The run loses orthogonality completely, \(\lVert V_k^\mathsf{T}V_k-I\rVert \approx1.1\). Nevertheless, the dilated residual satisfies \(\lVert \widehat{E}\rVert \approx2\cdot10^{-7}\,\lVert A\rVert \), and the completion of Lemma 8 produces an orthogonal \(\widehat{U}\) and a tridiagonal \(\widetilde{T}=\widehat{U}^\mathsf{T}(\widehat{A}+\vec\Delta)\widehat{U}\) whose leading \(k\times k\) block reproduces the computed \(T_k\) to working double precision, with \(\lVert \widehat{U}^\mathsf{T}\widehat{A}\widehat{U}-\widetilde{T}\rVert =\lVert \vec\Delta\rVert \approx3\cdot10^{-7}\). Every eigenvalue of \(\widetilde{T}\) therefore lies within \(\lVert \vec\Delta\rVert \) of \(\operatorname{spec}(A)\). Lanczos in exact arithmetic (double precision with full reorthogonalization) applied to \((\widetilde{T},e_1)\) reproduces every computed coefficient \(\alpha_j\) and \(\beta_j\), while applied to \((A,\bar{v}_1)\) it departs from the computed run at order \(\lVert A\rVert \) as orthogonality is lost.
Figure 3 compares \(\widetilde{T}\) with \(A\) through their starting vectors. By item (ii) of Theorem 5, the weights of \(e_1\) in an eigenbasis of \(\widehat{U}^\mathsf{T}\widehat{A}\widehat{U}\) agree exactly with those of \(\bar{v}_1\) in an eigenbasis of \(A\); the figure shows the corresponding comparison for \(\widetilde{T}\). Each eigenvalue of \(A\) is matched in \(\widetilde{T}\) by a cluster of \(k\) eigenvalues within \(\lVert \vec\Delta\rVert \), and the total weight of \(e_1\) on each cluster reproduces the weight of \(\bar{v}_1\) on the corresponding eigenvalue of \(A\), here to within \(4\cdot10^{-6}\).
The kernel and normalization model yields the perturbed recurrence, near-normalization, and local products. These feed both the Paige Ritz conclusions and the physical-dilation proof of the Greenbaum model. All dependence on implementation accuracy is expressed through \(\varepsilon\); all discarded factors are polynomial only in the iteration count. The backward construction uses only the local Paige package, while the global Paige theorems independently explain when and how orthogonality is lost and why every computed Ritz value remains meaningful.
AI tools (Claude Code with Fable 5 / Opus 5 and Codex with GPT-5.6 Sol) were used throughout the preparation of this note. The author assumes full responsibility for the content.
\printbibliography
[1] The Lanczos and Conjugate Gradient Algorithms: From Theory to Finite Precision Computations Society for Industrial and Applied Mathematics 2006 10.1137/1.9780898718140
[2] Iterative Methods for Solving Linear Systems Society for Industrial and Applied Mathematics 1997 10.1137/1.9781611970937
[3] Randomized algorithms for low-rank matrix approximation: Design, analysis, and applications 2023
[4] The Lanczos algorithm for matrix functions: a handbook for scientists 2024
[5] An iteration method for the solution of the eigenvalue problem of linear differential and integral operators Journal of research of the National Bureau of Standards 1950 45 255-282
[6] The Computation of Eigenvalues and Eigenvectors of Very Large Sparse Matrices University of London 1971
[7] Computational Variants of the Lanczos Method for the Eigenproblem IMA Journal of Applied Mathematics 1972 10 3 373–381 10.1093/imamat/10.3.373
[8] Error Analysis of the Lanczos Algorithm for Tridiagonalizing a Symmetric Matrix IMA Journal of Applied Mathematics 1976 18 3 341–349 10.1093/imamat/18.3.341
[9] Accuracy and effectiveness of the Lanczos algorithm for the symmetric eigenproblem Linear Algebra and its Applications 1980 34 235–258 Dec 10.1016/0024-3795(80)90167-6
[10] Behavior of slightly perturbed Lanczos and conjugate-gradient recurrences Linear Algebra and its Applications 1989 113 7–63 Feb 10.1016/0024-3795(89)90285-1
[11] Predicting the Behavior of Finite Precision Lanczos and Conjugate Gradient Computations SIAM Journal on Matrix Analysis and Applications 1992 13 1 121–137 Jan 10.1137/0613011
[12] The Lanczos and conjugate gradient algorithms in finite precision arithmetic Acta Numerica 2006 15 471–542 May 10.1017/s096249290626001x
[13] The Symmetric Eigenvalue Problem Society for Industrial and Applied Mathematics 1998 10.1137/1.9781611971163
[14] Accuracy and Stability of Numerical Algorithms Society for Industrial and Applied Mathematics 2002 Philadelphia Second
[15] Error Bounds in the Simple Lanczos Procedure for Computing Functions of Symmetric Matrices and Eigenvalues Comput. Math. Math. Phys. 1991 31 7 20–30 7
[16] The Simple Lanczos Procedure: Estimates of the Error of the Gauss Quadrature Formula and Their Applications Comput. Math. Math. Phys. 1996 36 11 1481–1492 1
[17] Stability of the Lanczos Method for Matrix Function Approximation 1605–1624 Society for Industrial and Applied Mathematics 2018 Jan 10.1137/1.9781611975031.105
[18] Accuracy of the Lanczos Process for the Eigenproblem and Solution of Equations SIAM Journal on Matrix Analysis and Applications 2019 40 4 1371–1398 Jan 10.1137/17m1133725
[19] On Stabilization and Convergence of Clustered Ritz Values in the Lanczos Method SIAM Journal on Matrix Analysis and Applications 2005 27 3 891–908 Jan 10.1137/040608908
[20] On the real convergence rate of the conjugate gradient method Linear Algebra and its Applications 1991 154–156 535–549 10.1016/0024-3795(91)90393-B