Skip to content

Instantly share code, notes, and snippets.

@ericphanson
Last active Jan 19, 2019
Embed
What would you like to do?
Simultaneous diagonalization for complex commuting normal matrices
# 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.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment