Simultaneous diagonalization for complex commuting normal matrices

simdiag.jl

# Adapted from
# https://uk.mathworks.com/matlabcentral/fileexchange/46794-simdiag-m
# which has the following license:
#    Copyright (c) 2009, Christian B. Mendl
#    All rights reserved.
#
#    Redistribution and use in source and binary forms, with or without
#    modification, are permitted provided that the following conditions are met:
#
#    * Redistributions of source code must retain the above copyright notice, this
#      list of conditions and the following disclaimer.
#
#    * Redistributions in binary form must reproduce the above copyright notice,
#      this list of conditions and the following disclaimer in the documentation
#      and/or other materials provided with the distribution
#    * Neither the name of  nor the names of its
#      contributors may be used to endorse or promote products derived from this
#      software without specific prior written permission.
#    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
#    AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
#    IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
#    DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
#    FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
#    DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
#    SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
#    CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
#    OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
#    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

using LinearAlgebra

# Exact minimizer of
# |s c conj(a0) - c^2 conj(a21) + s^2 conj(a12)|^2 + |s c a0 + s^2 a21 - c^2 a12|^2
# Refer to
#	H. H. Goldstine and L. P. Horwitz, A Procedure for the
#	Diagonalization of Normal Matrices, J. ACM [1959]
function calc_min(a0, a21, a12)
    u = real(a0)
    v = imag(a0)

    tmp = (a21 + conj(a12)) / 2
    r = abs(tmp); beta = angle(tmp)

    tmp = (a21 - conj(a12)) / 2
    s = abs(tmp); gamma = angle(tmp)

    nu = beta - gamma
    sin_nu = sin(nu)
    cos_nu = cos(nu)

    L = u * v - 4 * r * s * sin_nu
    M = u^2 - v^2 + 4 * (r^2 - s^2)

    A = L * (r^2 - s^2) * sin_nu + M * r * s
    B = L * (r^2 + s^2) * cos_nu
    C = L * (r^2 - s^2) + M * r * s * sin_nu

    tmp = r * s * cos_nu * sqrt(M^2 + 4 * L^2)
    phi = (atan(-A*C+B*tmp,B*C+A*tmp)-beta-gamma)/2;

    r_cos_ba = r * cos(beta + phi)
    s_sin_ga = s * sin(gamma + phi)
    kappa = u^2 + v^2 - 4 * (r_cos_ba^2 + s_sin_ga^2)
    lambda = 4 * (u * r_cos_ba + v * s_sin_ga)
    theta = -atan(-lambda, kappa) / 4

    c = cos(theta)
    s = exp( im * phi) * sin(theta)
    return ComplexF64(c), ComplexF64(s)
end

# Approximate minimizer of
# |s c conj(v[:,1]) - c^2 conj(v[:,2]) + s^2 conj(v[:,3])|^2 + |s c v[:,1] + s^2 v[:,2] - c^2 v[:,3]|^2
# for c = cos(theta), s = exp(i phi) sin(theta)
function approx_min(v)
    target(c, s, v) = norm(s * c * conj.(v[:,1]) - c^2 * conj.(v[:,2]) + s^2 * conj.(v[:,3]), 2)^2 + norm(s * c * v[:,1] + s^2 * v[:,2] - c^2 * v[:,3], 2)^2

    c, s = calc_min(v[1,1], v[1,2], v[1,3])
    m = target(c, s, v)
    for j = 2:size(v, 1)
        c1, s1 = calc_min(v[j,1], v[j,2], v[j,3])
        x = target(c1, s1, v)
        if x < m
            m = x
            c = c1; s = s1
        end
    end

    return ComplexF64(c), ComplexF64(s)
end

## A = A*R, R = R[j,k,c,s]
function timesR!(A, j, k, c, s)
# A = A*R
    A[:,[j,k]] = [(c * A[:,j] + s * A[:,k]) (-conj(s) * A[:,j] + conj(c) * A[:,k])]
end

## A = R'*A*R, R = R[j,k,c,s]
function rotate!(A, j, k, c, s)
        A[:,[j,k]] = [(c * A[:,j] + s * A[:,k]) (-conj(s) * A[:,j] + conj(c) * A[:,k])];	# A = A*R
        A[[j,k],:] = [(conj(c) * A[j,:] + conj(s) * A[k,:]) (-s * A[j,:] + c * A[k,:])]';	# A = R'*A
end


"""
    simdiag(list; max_iter::Int64 = 100)

Given `list`, a vector of pairwise commuting normal matrices,
returns `Q`, a unitary matrix containing the simultaneous eigenvectors,
and the list of matrices in their simultaneous eigenbasis.

Reference:
    Angelika Bunse-Gerstnert, Ralph Byers, and Volker Mehrmann
 	Numerical Methods for Simultaneous Diagonalization
 	SIAM J. Matrix Anal. Appl. Vol. 14, No. 4, pp. 927-949, October 1993
"""
function simdiag(list; max_iter::Int64 = 100)
    list = [ComplexF64.(A) for A in list]
    tol = eps()^1.5
    n = size(list[1],1)

    Q =  Matrix{ComplexF64}(I,n, n)

    calc_off2(A) = norm(A - diagm( 0 => diag(A)))^2

    num_iter = 0

    off2 = sum(calc_off2.(list))
    nscale = sum(norm.(list))

    while off2 > tol * nscale
        for j = 1:n
            for k = (j + 1):n
                v = zeros(ComplexF64, length(list), 3)
                for m = eachindex(list)
                    v[m,:] = [list[m][j, j] - list[m][k, k], list[m][j, k], list[m][k, j]]
                end
                c, s = approx_min(v)
                timesR!(Q, j, k, c, s)
                for m = eachindex(list)
                    rotate!(list[m], j, k, c, s)
                end
            end
        end

        off2 = sum(calc_off2.(list))

        num_iter += 1
        if num_iter > max_iter
            println("Exiting: Maximum number of iterations ($(max_iter)) exceeded. Current relative error: $(off2 / nscale).")
            return Q, list
        end
    end

    return Q, list
end



# Testing with the 513 code (for quantum error correction)
using Test

@testset "513-code test" begin

    ⊗ = kron
    X = ComplexF64[0 1; 1 0]
    Z = ComplexF64[1 0; 0 -1]
    I2 = Matrix{ComplexF64}(I, 2, 2)

    g = [   X ⊗ Z ⊗ Z ⊗ X ⊗ I2,
            I2 ⊗ X ⊗ Z ⊗ Z ⊗ X,
            X ⊗ I2 ⊗ X ⊗ Z ⊗ Z,
            Z ⊗ X ⊗ I2 ⊗ X ⊗ Z ]

    # Make sure they do simultaneously commute
    for gx in g, gy in g
     @assert gx*gy ≈ gy*gx
    end

    Q, gdiag = simdiag(g);

    for x in g
        y = Q' * x * Q
        @test y ≈ diagm( 0 => diag(y) )
    end
end
# test passes

function randunitary(d)
   rg1 = randn(d,d)
   rg2 = randn(d,d)
   RG = rg1 + im*rg2
   Q,R = qr(RG);
   r = diag(R)
   L = diagm(0 => r./abs.(r));
   return Q*L
end

function rand_commuting_normals(d, n)
    evals = [ rand(ComplexF64, d) for _ = 1:n ]
    U = randunitary(d)
    U, [U * diagm( 0 => ev ) * U' for ev in evals]
end



@testset "Random normal matrices test" begin
    d = 5
    n = 4
    U, normals = rand_commuting_normals(d, n)

    # Make sure they do simultaneously commute
    for x in normals, y in normals
     @assert x*y ≈ y*x
    end

    Q, list = simdiag(normals; max_iter = 2000);

    for x in list
        y = Q' * x * Q
        @test y ≈ diagm( 0 => diag(y) )
    end
end

# Test fails...
# Increasing max_iter does not seem to help.