Remarks on the Unitary Triangularization of a Nonsymmetric Matrix

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 ‖ <italic>u</italic> ‖ = √2, then the matrix <italic>I</italic> - <italic>uu</italic><supscrpt>*</supscrpt> is unitary, as may be readily verified. The following lemma is a modification of the one give in [1]. LEMMA. <italic>Let a</italic> ≠ 0 <italic>be an arbitrary vector and let v be an arbitrary unit vector. Then there exists a vector u with</italic> ‖ <italic>u</italic> ‖ = √2 <italic>and a scalar &zgr; with</italic> | <italic>&zgr;</italic> | = 1 <italic>such that</italic> (<italic>I</italic> - <italic>uu</italic><supscrpt>*</supscrpt>)<italic>a</italic> = <italic>&zgr;</italic> ‖ <italic>a</italic> ‖ <italic>v</italic>. (1) PROOF: Letting <italic>α</italic> = ‖ <italic>a</italic> ‖ and <italic>μ</italic> = <italic>u</italic><supscrpt>*</supscrpt><italic>a</italic>, (1) may be written <italic>a</italic> - <italic>α&zgr;v</italic> = <italic>μu</italic>. (2) Multiply by <italic>a</italic><supscrpt>*</supscrpt> gives <italic>α</italic><supscrpt>2</supscrpt> - <italic>α&zgr;a</italic><supscrpt>*</supscrpt><italic>v</italic> = <italic>μa</italic><supscrpt>*</supscrpt><italic>u</italic> = ‖ <italic>μ</italic> ‖<supscrpt>2</supscrpt>. (3) It follows that <italic>&zgr;a</italic><supscrpt>*</supscrpt><italic>v</italic> is real. Assuming for the moment that <italic>a</italic><supscrpt>*</supscrpt><italic>v</italic> ≠ 0, we write it in polar form <italic>a</italic><supscrpt>*</supscrpt><italic>v</italic> = <italic>rw</italic>, <italic>r</italic> > 0, ‖ <italic>w</italic> ‖ = 1. Then the fact that <italic>&zgr;rw</italic> is real implies that <italic>&zgr;</italic> = ± <italic>w</italic>. (4) Substituting into (3) gives ‖ <italic>μ</italic> ‖<supscrpt>2</supscrpt> = <italic>α</italic><supscrpt>2</supscrpt> ∓ <italic>αr</italic>. (5) We now set, arbitrarily, arg (<italic>μ</italic>) = 0. Then <italic>μ</italic> = √<italic>α</italic>(<italic>α</italic> ∓ <italic>r</italic>). (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 <italic>u</italic> and <italic>&zgr;</italic>: <italic>α</italic> = ‖ <italic>a</italic> ‖ (7) <italic>r</italic> = ‖ <italic>a</italic><supscrpt>*</supscrpt><italic>v</italic> ‖ (8) <italic>&zgr;</italic> = -<italic>a</italic><supscrpt>*</supscrpt><italic>v</italic>/<italic>r</italic> (9) <italic>μ</italic> = √ <italic>α</italic>(<italic>α</italic> + <italic>r</italic>) (10) <italic>u</italic> = 1/<italic>μ</italic>(<italic>a</italic> - <italic>&zgr;αv</italic>). (11) The case <italic>a</italic><supscrpt>*</supscrpt><italic>v</italic> = 0 may be handled if, instead of using (9), we let <italic>&zgr;</italic> be an arbitrary number with ‖ <italic>&zgr;</italic> ‖ = 1. It is easily verified that the formulas thus modified will still work. The computation requires 3 square roots to compare <italic>α</italic>, <italic>r</italic>, and <italic>μ</italic>. A slight modification [1] permits one to avoid the root required to compute <italic>μ</italic>. In the real case, no root is required to compute <italic>r</italic>. 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 <italic>r</italic> = <italic>a</italic><supscrpt>*</supscrpt><italic>v</italic>, and (10) by <italic>μ</italic> = √<italic>α</italic>(<italic>α</italic> - <italic>r</italic>). If <italic>μ</italic> is computed this way and if <italic>r</italic> is positive and near <italic>α</italic> (as is the case when <italic>v</italic> is near <italic>a</italic>/<italic>α</italic>), cancellation of significant digits will occur. This difficulty, and the need for making a special case when <italic>v</italic> is exactly <italic>a</italic>/<italic>α</italic>, is avoided in the present set of formulas.