Remarks on the Unitary Triangularization of a Nonsymmetric Matrix
D. D. Morrison · Journal of the ACM · 1960
In [1], A. Householder described a method for the unitary triangularization of a matrix. The formulas given there are valid for the real case. In this note we describe the modifications to handle the complex case and also point out a small modification in the real case which will improve the numerical accuracy of the method. At first we are concerned with a complex vector space. The basic tool is the fact that if ‖ u ‖ = √2, then the matrix I - uu * is unitary, as may be readily verified. The following lemma is a modification of the one give in [1]. LEMMA. Let a ≠ 0 be an arbitrary vector and let v be an arbitrary unit vector. Then there exists a vector u with ‖ u ‖ = √2 and a scalar ζ with | ζ | = 1 such that ( I - uu * ) a = ζ ‖ a ‖ v . (1) PROOF: Letting α = ‖ a ‖ and μ = u * a , (1) may be written a - αζv = μu . (2) Multiply by a * gives α 2 - αζa * v = μa * u = ‖ μ ‖ 2 . (3) It follows that ζa * v is real. Assuming for the moment that a * v ≠ 0, we write it in polar form a * v = rw , r > 0, ‖ w ‖ = 1. Then the fact that ζrw is real implies that ζ = ± w . (4) Substituting into (3) gives ‖ μ ‖ 2 = α 2 ∓ αr . (5) We now set, arbitrarily, arg ( μ ) = 0. Then μ = √ α ( α ∓ r ). (6) Next, we select the negative sign in (4) in order to avoid the subtraction of two positive quantities in (6), since such a subtraction may give rise to numerical difficulties. Collecting the formulas, we see that the following sequence of computations will produce the required u and ζ : α = ‖ a ‖ (7) r = ‖ a * v ‖ (8) ζ = - a * v / r (9) μ = √ α ( α + r ) (10) u = 1/ μ ( a - ζαv ). (11) The case a * v = 0 may be handled if, instead of using (9), we let ζ be an arbitrary number with ‖ ζ ‖ = 1. It is easily verified that the formulas thus modified will still work. The computation requires 3 square roots to compare α , r , and μ . A slight modification [1] permits one to avoid the root required to compute μ . In the real case, no root is required to compute r . Now consider the case of a real vector space. The formulas given [1] for this case are essentially the same as ours except that (8) is replaced by r = a * v , and (10) by μ = √ α ( α - r ). If μ is computed this way and if r is positive and near α (as is the case when v is near a / α ), cancellation of significant digits will occur. This difficulty, and the need for making a special case when v is exactly a / α , is avoided in the present set of formulas.