J -Orthogonal Matrices: Properties and
Generation
Higham, Nicholas J.
2003
MIMS EPrint: 2006.69
Manchester Institute for Mathematical Sciences
School of Mathematics
The University of Manchester
Reports available from: http://eprints.maths.manchester.ac.uk/
And by contacting: The MIMS Secretary
School of Mathematics
The University of Manchester
Manchester, M13 9PL, UK
ISSN 1749-9097
J-ORTHOGONAL MATRICES: PROPERTIES AND GENERATION
∗
NICHOLAS J. HIGHAM
†
Abstract. A real, square matrix Q is J-orthogonal if Q
T
JQ = J, where the signature matrix
J = diag(±1). J-orthogonal matrices arise in the analysis and numerical solution of various matrix
problems involving indefinite inner products, including, in particular, the downdating of Cholesky
factorizations. We present techniques and tools useful in the analysis, application and construction
of these matrices, giving a self-contained treatment that provides new insights. First, we define and
explore the properties of the exchange operator, which maps J-orthogonal matrices to orthogonal
matrices and vice versa. Then we show how the exchange operator can be used to obtain a hyperbolic
CS decomposition of a J-orthogonal matrix directly from the usual CS decomposition of an orthogonal
matrix. We employ the decomposition to derive an algorithm for constructing random J-orthogonal
matrices with specified norm and condition numb er. We also give a short proof of the fact that
J-orthogonal matrices are optimally scaled under two-sided diagonal scalings. We introduce the
indefinite polar decomposition and investigate two iterations for computing the J-orthogonal polar
factor: a Newton iteration involving only matrix inversion and a Schulz iteration involving only
matrix multiplication. We show that these iterations can be used to J-orthogonalize a matrix that
is not too far from being J-orthogonal.
Key words. J-orthogonal matrix, exchange operator, gyration operator, sweep operator, princi-
pal pivot transform, hyperbolic CS decomposition, two-sided scaling, indefinite least squares problem,
hyperbolic QR factorization, indefinite polar decomposition, Newton’s method, Schulz iteration
AMS subject classifications. 65F30, 15A18
1. Introduction. A matrix Q ∈ R
n×n
is J-orthogonal if
Q
T
JQ = J,(1.1)
where J = diag(±1) is a signature matrix. Clearly, Q is nonsingular and QJQ
T
= J.
This type of matrix arises in hyperbolic problems, that is, problems where there is
an underlying indefinite inner product or weight matrix. We give two examples to
illustrate the utility of J-orthogonal matrices.
First consider the downdating problem of computing the Cholesky factorization of
a positive definite matrix C = A
T
A − B
T
B, where A ∈ R
p×n
(p ≥ n) and B ∈ R
q ×n
.
This task arises when solving a regression problem after some of the rows (namely
those of B) of the data matrix are removed, and A in this case is usually upper
triangular. Numerical stability considerations dictate that we should avoid explicit
formulation of C. If we can find a J-orthogonal matrix Q such that
Q
A
B
=
R
0
,(1.2)
with J = diag(I
p
, −I
q
) and R ∈ R
n×n
upper triangular, then
C =
A
B
T
J
A
B
=
A
B
T
Q
T
JQ
A
B
= R
T
R,
∗
Numerical Analysis Report 408, Manchester Centre for Computational M athematics, September
2002. Revised March 2003.
†
Department of Mathematics, University of Manchester, Manchester, M13 9PL, England
(higham@ma.man.ac.uk, http://www.ma.man.ac.uk/~higham/). This work was supported by Engi-
neering and Physical Sciences Research Council grant GR/R22612.
1
2 NICHOLAS J. HIGHAM
so R is the desired Cholesky factor. The factorization (1.2) is a hyperbolic QR fac-
torization; for details of how to compute it see, for example, [1].
A second example where J-orthogonal matrices play a key role is in the solution
of the symmetric definite generalized eigenproblem Ax = λBx, where A and B are
symmetric, some linear combination of them is positive definite, and B is nonsingular.
Through the use of a congruence transformation (for example by using a block LDL
T
decomposition of B followed by a diagonalization of the block diagonal factor [38]) the
problem can be reduced to
e
Ax = λJx, for some signature matrix J = diag (±1). If we
can find a J-orthogonal Q such that Q
T
e
AQ = D = diag (d
i
) then the eigenvalues are
the diagonal elements of JD; such a Q can be constructed using a Jacobi algorithm
of Veseli´c [40].
In addition to these practical applications, J-orthogonal matrices are of significant
theoretical interest. For example, they play a fundamental role in the study of J-
contractive matrices [30], which are matrices X for which XJX
T
≤ J, where A ≥ 0
denotes that the symmetric matrix A is p ositive semidefinite.
A matrix Q ∈ R
n×n
is (J
1
, J
2
)-orthogonal if
Q
T
J
1
Q = J
2
,(1.3)
where J
1
= diag (±1) and J
2
= diag (±1) are signature matrices having the same
inertia. (J
1
, J
2
)-orthogonal matrices are also known as hyperexchange matrices and
J-orthogonal matrices as hypernormal matrices [2]. Since J
1
and J
2
in (1.3) have the
same inertia, J
2
= P J
1
P
T
for some permutation matrix P , and hence (QP )
T
J
1
(QP ) =
J
1
. A (J
1
, J
2
)-orthogonal matrix is therefore simply a column permutation of a J
1
-
orthogonal matrix, and so for the purposes of this work we can restrict our attention
to J-orthogonal matrices. An application in which (J
1
, J
2
)-orthogonal matrices arise
with J
1
and J
2
generally different is the HR algorithm of Brebner and Grad [5] and
Bunse-Gerstner [6] for solving the standard eigenvalue problem for J-symmetric ma-
trices. A matrix A ∈ R
n×n
is J-symmetric if AJ is symmetric, or, equivalently, if
JA, JA
T
or A
T
J is symmetric. Given a J
0
-symmetric matrix A, the kth stage of the
unshifted HR algorithm consists of factoring A
k
= H
k
R
k
, where H
k
is (J
k
, J
k+1
)-
orthogonal with J
k+1
= H
T
k
J
k
H
k
and R
k
is upper triangular, and then setting
A
k+1
= R
k
H
k
. Computational details and convergence properties of the algorithm
can be found in [6].
Unlike the subclass of orthogonal matrices, J-orthogonal matrices can be arbi-
trarily ill conditioned. This poses interesting questions and difficulties in the design,
analysis and testing of algorithms and motivates our attempt to gain a better under-
standing of the class of J-orthogonal matrices.
The purpose of this paper is threefold. First we collect some interesting and not
so well-known properties of J-orthogonal matrices. In particular, we give a new proof
of the hyperbolic CS decomposition via the usual CS decomposition by exploiting the
exchange operator. The exchange operator is a tool that has found use in several areas
of mathematics and is known by several different names; we give a brief survey of its
properties and its history. We also give a new proof of the fact that J-orthogonal
matrices are optimally scaled under two-sided diagonal scalings. Our second aim is to
show how to generate random J-orthogonal matrices with specified singular values,
and in particular with specified norms and condition numbers—a capability that is
very useful for constructing test data for problems with an indefinite flavour. Finally,
we investigate two Newton iterations for computing a J-orthogonal matrix, one involv-
ing only matrix inversion, the other only matrix multiplication. Both iterations are
J-ORTHOGONAL MATRICES 3
shown to conver ge to the J-orthogonal factor in a certain indefinite polar decomposi-
tion under suitable conditions. Analogously to the case of orthogonal matrices and the
corresponding Newton iterations [15], [20], we show that these Newton iterations can
be used to J-orthogonalize a matrix that is not too far from being J-orthogonal. An
application is to the situation where a matrix that should be J-orthogonal turns out
not to be because of rounding or other errors and it is desired to J-orthogonalize it.
J-orthogonal matrices, and hyperbolic problems in general, are the subject of
much recent and current research, covering both theory and algorithms. This paper
provides a self-contained treatment that highlights some techniques and tools useful
in the analysis and application of these matrices; the treatment should also be of more
general interest.
Throughout, we take J to have the form
J =
I
p
0
0 −I
q
, p + q = n,(1.4)
and we use exclusively the 2-norm: kAk
2
= max
x6=0
kAxk
2
/kxk
2
, where kxk
2
2
= x
T
x.
2. The exchange operator. Let A ∈ R
n×n
and consider the system
y =
1
p y
1
q y
2
=
p q
p A
11
A
12
q A
21
A
22
1
x
1
x
2
p
q
= Ax,(2.1)
where A
11
is nonsingular. We use this partitioning of A throughout the section. By
solving the first equation in (2.1) for x
1
and then eliminating x
1
from the second
equation we obtain
x
1
y
2
= exc(A)
y
1
x
2
,(2.2)
where
exc(A) =
A
−1
11
−A
−1
11
A
12
A
21
A
−1
11
A
22
− A
21
A
−1
11
A
12
.
We call exc the exchange operator, since it exchanges x
1
and y
1
in (2.1). Note that
the (2,2)-block of exc(A) is the Schur complement of A
11
in A. The definition of
the exchange operator can be generalized to allow the “pivot matrix” A
11
to be any
principal submatrix, but for our purposes this extra level of generality is not necessary.
It is easy to see that the exchange operator is involutary,
exc(exc(A)) = A,(2.3)
and moreover that
exc(JAJ) = Jexc(A)J = exc(A
T
)
T
.(2.4)
This last identity shows that J is naturally associated with exc.
We first address the nonsingularity of exc(A). The blo ck LU factorization
exc(A) =
I 0
A
21
A
22
A
−1
11
−A
−1
11
A
12
0 I
=
I 0
A
21
A
22
A
11
A
12
0 I
−1
≡ LR
−1
(2.5)
4 NICHOLAS J. HIGHAM
will be useful.
Lemma 2.1. Let A ∈ R
n×n
with A
11
nonsingular. Then exc(A) is nonsingular if
and only if A
22
is nonsingular. If A is nonsingular and exc(A
−1
) exists then exc(A)
is nonsingular and
exc(A)
−1
= exc(A
−1
).(2.6)
Proof. For A ∈ R
n×n
with A
11
nonsingular, the block LU factorization (2.5)
makes clear that exc(A) is nonsingular if and only if A
22
is nonsingular. The last part
is obtained by rewriting (2.1) as x = A
−1
y and deriving the corresponding analogue
of (2.2):
y
1
x
2
= exc(A
−1
)
x
1
y
2
.(2.7)
It follows from (2.7) that for any x
1
and y
2
there is a unique x
2
and y
1
, which implies
from (2.2) that exc(A) is nonsingular and exc(A)
−1
= exc(A
−1
).
Note that either of A and exc(A) can be singular without the other being singular,
as shown by the examples with p = q = 1,
A =
1 1
1 0
, exc(A) =
1 −1
1 −1
, A =
1 1
1 1
, exc(A) =
1 −1
1 0
.
For completeness, we mention that for both exc(A) and exc(A
−1
) to exist and be
nonsingular, it is necessary and sufficient that A, A
11
and A
22
be nonsingular.
The reason for our interest in the exchange operator is that it maps J-orthogonal
matrices to orthogonal matrices and vice versa. Note that J-orthogonality of A implies
that A
T
11
A
11
= I + A
T
21
A
21
and hence that A
11
is nonsingular and exc(A) exists, but
if A is orthogonal A
11
can be singular.
Theorem 2.2. Let A ∈ R
n×n
. If A is J-orthogonal then exc(A) is orthogonal.
If A is orthogonal and A
11
is nonsingular then exc(A) is J-orthogonal.
Proof. Proving the result by working directly with exc(A) involves some laborious
algebra. A more elegant proof involving quadratic forms is given by Stewart and
Stewart [36, sec. 2]. We give another proof, suggested by Chris Paige. Assume first
that A is orthogonal with A
11
nonsingular. Then exc(A
T
) = exc(A
−1
) exists and
Lemma 2.1 shows that exc(A) is nonsingular and exc(A)
−1
= exc(A
−1
) = exc(A
T
).
Hence, using (2.4),
I = exc(A
T
)exc(A) = Jexc(A)
T
J · exc(A),
which shows that exc(A) is J-orthogonal.
If A is J-orthogonal then, as noted above, A
11
is nonsingular. Also JA
T
J = A
−1
and so from Lemma 2.1, exc(JA
T
J) = exc(A
−1
) = exc(A)
−1
. But (2.4) shows that
exc(JA
T
J) = exc(A)
T
, and we conclude that exc(A) is orthogonal.
As an example of a result of a different flavour, we give the following generalization
to arbitrary p of a result obtained by Duffin, Hazony, and Morrison [10] for p = 1.
Theorem 2.3. Let A ∈ R
n×n
with A
11
nonsingular. Then exc(A) + exc(A)
T
is
congruent to A + A
T
.
Proof. Using (2.5) we have
exc(A) + exc(A)
T
= LR
−1
+ R
−T
L
T
= R
−T
(R
T
L + L
T
R)R
−1
.