How many arithmetic operations does matlab require to determine the Schur decomposition?

Consider a matrix A of size n.

I would like to determine the square root of this matrix which can be done like this:

n = 10; % variable size 
A = rand(n);  % random square matrix A of size n by n
A_sqrt = sqrtm(A); % square root of matrix A

Inside the sqrtm command, matlab requires the Schur decomposition for determining the square root of a matrix.

[S, T] = schur(A); % with A = S*T*S and S*S' = I and T = (quassi) upper triangular matrix 

To determine the speed of sqrtm, I would like an expression of the amount of required distinguisable operations (summation, subtraction, multiplication, division and square root) expressed in the matrix size n. To get this expression, I would like to know how Matlab determines the Schur decomposition of a matrix.

I read that the first step is to determine the upper Hessenberg form H by means of

[G,H] = hess(A); % with A = G*H*G' and G*G' = I

After this, a QR decomposition of H is executed

[Q,R] = qr(H); % with H = Q*R and Q*Q' = I and R = upper triangular matrix

How does one continue from here to a Schur decomposition?

If this is not the way how Matlab determines a Schur decomposition, what is it?

Answer 1

[Christine Tobler]

The second step in computing the Schur decomposition is the QR algorithm, not the QR decomposition (unfortunately the names are very similar). The Hessenberg form is computed in advance to allow the QR algorithm to be applied implicitly. There are other ways of computing the Schur decomposition, but the QR algorithm is the most standard and simplest.

It's not very practical to count the number of operations in this algorithm, particularly because the QR algorithm is iterative, so the number of passes through the data depends on the convergence speed, which again depends on the eigenvalues that are used. In practice, the locality of the code (how many times a block of memory is accessed) is also more relevant to performance than simply counting how many times each operation is done.

In practice, I would normally say that Schur has complexity O(n^3), and check how the performance by getting time measurements on a specific machine.

[Jip Reinders]

Thanks Christine, this clarifies my problem.

I do have a follow up question: Does the problem changes when matrix A is symmetric? the matrix A is defined as

N = 4;
n = 2; % amount of states
m = 2; % amount of measurements (n = m)
X = rand(n,N);
C = rand(m,n); % Observer matrix in state space model
V = rand(m,m); V = V + V’; % symmetric measurement noise covariance matrix 
I = eye(N);
A = inv(I + X’*C’*V*C*X); % this should result in a symmetric matrix

in that case it would also be possible to use

[E, D] = eig(A) % with V*V’ = I and D = diagonal matrix

the reason behind the requirement of the decomposition is because I need to determine the square root of matrix A sqrtm(A) which can be done as

R = sqrt(diag(D));
A_sqrt  = (E.*R.’)*E‘
A_sqrt = (A_sqrt + A_sqrt’)/2 % I took this part from the code behind sqrtm in Matlab

second follow up question: Does the problem changes when matrix A is symmetric possitive definite?

thanks in advance for responding

[Christine Tobler]

The algorithm used changes when the matrix A is symmetric: The output of hess is a tridiagonal matrix, which makes the implicit QR algorithm much cheaper to apply. Also, this algorithm will guarantee that the resulting eigenvalues and eigenvectors are real, which is not the case for non-symmetric inputs.

Note that EIG will check if the input matrix is exactly symmetric (issymmetric(A)) and only use the symmetric algorithm in that case. The formula inv(I + X’*C’*V*C*X) will probably result in a nonsymmetric matrix, and passing it into EIG as (A + A')/2 should speed up that computation.