# Neat Randomized Algorithms: RandDiag for Rapidly Diagonalizing Normal Matrices

Consider two complex-valued square matrices and . The first matrix is Hermitian, being equal to its conjugate transpose . The other matrix is non-Hermitian, . Let’s see how long it takes to compute their eigenvalue decompositions in MATLAB:

>> A = randn(1e3) + 1i*randn(1e3); A = (A+A')/2;
>> tic; [V_A,D_A] = eig(A); toc % Hermitian matrix
Elapsed time is 0.415145 seconds.
>> B = randn(1e3) + 1i*randn(1e3);
>> tic; [V_B,D_B] = eig(B); toc % non-Hermitian matrix
Elapsed time is 1.668246 seconds.

We see that it takes longer to compute the eigenvalues of the non-Hermitian matrix as compared to the Hermitian matrix . Moreover, the matrix of eigenvectors for a Hermitian matrix is a unitary matrix, .

There are another class of matrices with nice eigenvalue decompositions, normal matrices. A square, complex-valued matrix is normal if . The matrix of eigenvectors for a normal matrix is also unitary, . An important class of normal matrices are unitary matrices themselves. A unitary matrix is always normal since it satisfies .

Let’s see how long it takes MATLAB to compute the eigenvalue decomposition of a unitary (and thus normal) matrix:

>> U = V_A;                     % unitary, and thus normal, matrix
>> tic; [V_U,D_U] = eig(U); toc % normal matrix
Elapsed time is 2.201017 seconds.

Even longer than it took to compute an eigenvalue decomposition of the non-normal matrix ! Can we make the normal eigenvalue decomposition closer to the speed of the Hermitian eigenvalue decomposition?

Here is the start of an idea. Every square matrix has a Cartesian decomposition:

We have written as a combination of its Hermitian part and times its skew-Hermitian part . Both and are Hermitian matrices. The Cartesian decomposition of a square matrix is analogous to the decomposition of a complex number into its real and imaginary parts.

For a normal matrix , the Hermitian and skew-Hermitian parts commute, . We know from matrix theory that commuting Hermitian matrices are simultaneously diagonalizable, i.e., there exists such that and for diagonal matrices and . Thus, given access to such , has eigenvalue decomposition

Here’s a first attempt to turn this insight into an algorithm. First, compute the Hermitian part of , diagonalize , and then see if diagonalizes . Let’s test this out on a example:

>> C = orth(randn(2) + 1i*randn(2)); % unitary matrix
>> H = (C+C')/2;                     % Hermitian part
>> [Q,~] = eig(H);
>> Q'*C*Q                            % check to see if diagonal
ans =
-0.9933 + 0.1152i  -0.0000 + 0.0000i
0.0000 + 0.0000i  -0.3175 - 0.9483i

Yay! We’ve succeeded at diagonalizing the matrix using only a Hermitian eigenvalue decomposition. But we should be careful about declaring victory too early. Here’s a bad example:

>> C = [1 1i;1i 1]; % normal matrix
>> H = (C+C')/2;
>> [Q,~] = eig(H);
>> Q'*C*Q           % oh no! not diagonal
ans =
1.0000 + 0.0000i   0.0000 + 1.0000i
0.0000 + 1.0000i   1.0000 + 0.0000i

What’s going on here? The issue is that the Hermitian part for this matrix has a repeated eigenvalue. Thus, has multiple different valid matrices of eigenvectors. (In this specific case, every unitary matrix diagonalizes .) By looking at alone, we don’t know which matrix to pick which also diagonalizes .

He and Kressner developed a beautifully simple randomized algorithm called RandDiag to circumvent this failure mode. The idea is straightforward:

1. Form a random linear combination of the Hermitian and skew-Hermitian parts of , with standard normal random coefficients and .
2. Compute that diagonalizes .

That’s it!

To get a sense of why He and Kressner’s algorithm works, suppose that has some repeated eigenvalues and has all distinct eigenvalues. Given this setup, it seems likely that a random linear combination of and will also have all distinct eigenvalues. (It would take a very special circumstances for a random linear combination to yield two eigenvalues that are exactly the same!) Indeed, this intuition is a fact: With 100% probability, diagonalizing a Gaussian random linear combination of simultaneously diagonalizable matrices and also diagonalizes and individually.

MATLAB code for RandDiag is as follows:

function Q = rand_diag(C)
H = (C+C')/2; S = (C-C')/2i;
M = randn*H + randn*S;
[Q,~] = eig(M);
end

When applied to our hard example from before, RandDiag succeeds at giving a matrix that diagonalizes :

>> Q = rand_diag(C);
>> Q'*C*Q
ans =
1.0000 - 1.0000i  -0.0000 + 0.0000i
-0.0000 - 0.0000i   1.0000 + 1.0000i

For computing the matrix of eigenvectors for a unitary matrix, RandDiag takes 0.4 seconds, just as fast as the Hermitian eigendecomposition did.

>> tic; V_U = rand_diag(U); toc
Elapsed time is 0.437309 seconds.

He and Kressner’s algorithm is delightful. Ultimately, it uses randomness in only a small way. For most coefficients , a matrix diagonalizing will also diagonalize . But, for any specific choice of , there is a possibility of failure. To avoid this possibility, we can just pick and at random. It’s really as simple as that.

References: RandDiag was proposed in A simple, randomized algorithm for diagonalizing normal matrices by He and Kressner (2024), building on their earlier work in Randomized Joint Diagonalization of Symmetric Matrices (2022) which considers the general case of using random linear combinations to (approximately) simultaneous diagonalize (nearly) commuting matrices. RandDiag is an example of a linear algebraic algorithm that uses randomness to put the input into “general position”; see Randomized matrix computations: Themes and variations by Kireeva and Tropp (2024) for a discussion of this, and other, ways of using randomness to design matrix algorithms.

# Neat Randomized Algorithms: Randomized Cholesky QR

As a research area, randomized numerical linear algebra (RNLA) is as hot as ever. To celebrate the exciting work in this space, I’m starting a new series on my blog where I celebrate cool recent algorithms in the area. At some future point, I might talk about my own work in this series, but for now I’m hoping to use this series to highlight some of the awesome work being done by my colleagues.

Given a tall matrix with , its (economy-size) QR factorization is a decomposition of the form , where is a matrix with orthonormal columns and is upper triangular. QR factorizations are used to solve least-squares problems and as a computational procedure to orthonormalize the columns of a matrix.

Here’s an example in MATLAB, where we use QR factorization to orthonormalize the columns of a test matrix. It takes about 2.5 seconds to run.

>> A = randn(1e6, 1e2) * randn(1e2) * randn(1e2); % test matrix
>> tic; [Q,R] = qr(A,"econ"); toc
Elapsed time is 2.647317 seconds.

The classical algorithm for computing a QR factorization uses Householder reflectors and is exceptionally numerically stable. Since has orthonormal columns, is the identity matrix. Indeed, this relation holds up to a tiny error for the computed by Householder QR:

>> norm(Q'*Q - eye(1e2)) % || Q^T Q - I ||
ans =
7.0396e-14

The relative error is also small:

>> norm(A - Q*R) / norm(A)
ans =
4.8981e-14

Here is an alternate procedure for computing a QR factorization, known as Cholesky QR:

function [Q,R] = cholesky_qr(A)
R = chol(A'*A);
Q = A / R;       % Q = A * R^{-1}
end

This algorithm works by forming , computing its (upper triangular) Cholesky decomposition , and setting . Cholesky QR is very fast, about faster than Householder QR for this example:

>> tic; [Q,R] = cholesky_qr(A); toc
Elapsed time is 0.536694 seconds.

Unfortunately, Cholesky QR is much less accurate and numerically stable than Householder QR. Here, for instance, is the value of , about ten million times larger than for Householder QR!:

>> norm(Q'*Q - eye(1e2))
ans =
7.5929e-07

What’s going on? As we’ve discussed before on this blog, forming is typically problematic in linear algebraic computations. The “badness” of a matrix is measured by its condition number, defined to be the ratio of its largest and smallest singular values . The condition number of is the square of the condition number of , , which is at the root of Cholesky QR’s loss of accuracy. Thus, Cholesky QR is only appropriate for matrices that are well-conditioned, having a small condition number , say .

The idea of randomized Cholesky QR is to use randomness to precondition , producing a matrix that is well-conditioned. Then, since is well-conditioned, we can apply ordinary Cholesky QR to it without issue. Here are the steps:

1. Draw a sketching matrix of size ; see these posts of mine for an introduction to sketching.
2. Form the sketch . This step compresses the very tall matrix to the much shorter matrix of size .
3. Compute a QR factorization using Householder QR. Since the matrix is small, this factorization will be quick to compute.
4. Form the preconditioned matrix .
5. Apply Cholesky QR to to compute .
6. Set . Observe that , as desired.

MATLAB code for randomized Cholesky QR is provided below:1Code for the sparsesign subroutine can be found here.

function [Q,R] = rand_cholesky_qr(A)
S = sparsesign(2*size(A,2),size(A,1),8); % sparse sign embedding
R1 = qr(S*A,"econ"); % sketch and (Householder) QR factorize
B = A / R1; % B = A * R_1^{-1}
[Q,R2] = cholesky_qr(B);
R = R2*R1;
end

Randomized Cholesky QR is still faster than ordinary Householder QR, about faster in our experiment:

>> tic; [Q,R] = rand_cholesky_qr(A); toc
Elapsed time is 0.920787 seconds.

Randomized Cholesky QR greatly improves on ordinary Cholesky QR in terms of accuracy and numerical stability. In fact, the size of is even smaller than for Householder QR!

>> norm(Q'*Q - eye(1e2))
ans =
1.0926e-14

The relative error is small, too! Even smaller than for Householder QR in fact:

>> norm(A - Q*R) / norm(A)
ans =
4.0007e-16

Like many great ideas, randomized Cholesky QR was developed independently by a number of research groups. A version of this algorithm was first introduced in 2021 by Fan, Guo, and Lin. Similar algorithms were investigated in 2022 and 2023 by Balabanov, Higgins et al., and Melnichenko et al. Check out Melnichenko et al.‘s paper in particular, which shows very impressive results for using randomized Cholesky QR to compute column pivoted QR factorizations.

References: Primary references are A Novel Randomized XR-Based Preconditioned CholeskyQR Algorithm by Fan, Guo, and Lin (2021); Randomized Cholesky QR factorizations by Balabanov (2022); Analysis of Randomized Householder-Cholesky QR Factorization with Multisketching by Higgins et al. (2023); CholeskyQR with Randomization and Pivoting for Tall Matrices (CQRRPT) by Melnichenko et al. (2023). The idea of using sketching to precondition tall matrices originates in the paper A fast randomized algorithm for overdetermined linear least-squares regression by Rokhlin and Tygert (2008).