Simultaneous diagonalization (#37)

Implement simultaneous diagonalization of commuting operators.

src/QuantumOptics.jl

@@ -23,7 +23,7 @@ export bases, Basis, GenericBasis, CompositeBasis,
         nparticlebasis, BosonicNParticleBasis, FermionicNParticleBasis, nparticleoperator_1, nparticleoperator_2,
         metrics, tracedistance, tracedistance_general, entropy_vn,
         spectralanalysis, operatorspectrum, operatorspectrum_hermitian,
-                eigenstates, eigenstates_hermitian, groundstate,
+                eigenstates, eigenstates_hermitian, groundstate, simdiag,
         timeevolution_simple,
         timeevolution,
         cumulantexpansion,

src/spectralanalysis.jl

@@ -2,7 +2,7 @@ module spectralanalysis
 
 using ..states, ..operators, ..operators_dense, ..operators_sparse
 
-export operatorspectrum, operatorspectrum_hermitian, eigenstates, eigenstates_hermitian, groundstate
+export operatorspectrum, operatorspectrum_hermitian, eigenstates, eigenstates_hermitian, groundstate, simdiag
 
 
 """
@@ -158,4 +158,54 @@ Nmax (optional)
 """
 groundstate(H::Union{DenseOperator, SparseOperator}) = eigenstates_hermitian(H; Nmax=1)[1]
 
+
+"""
+Simultaneously diagonalize commuting Hermitian operators.
+
+This is done by diagonalizing the sum of the operators.
+The eigenvalues are computed by :math:`a = \\langle \\psi |A|\\psi\\rangle` and
+it is checked whether the eigenvectors fulfill the equation
+:math:`A|\\psi\\rangle = a|\\psi\\rangle`.
+
+Arguments
+---------
+
+Ops
+  Vector of operators (sparse or dense).
+
+Keyword arguments
+-----------------
+
+atol (optional)
+  kwarg of Base.isapprox specifying the tolerance of the approximate check
+  Default is 1e-14.
+
+rtol (optional)
+  kwarg of Base.isapprox specifying the tolerance of the approximate check
+  Default is 1e-14.
+"""
+function simdiag{T <: DenseOperator}(Ops::Vector{T}; atol::Real=1e-14, rtol::Real=1e-14)
+
+  # Check input
+  for A=Ops
+    if !isapprox(A.data, dagger(A).data; atol=atol, rtol=rtol)
+      error("Non-hermitian operator given!")
+    end
+  end
+
+  d, v = eig(sum(Ops).data)
+
+  evals = [Vector{Complex128}(length(d)) for i=1:length(Ops)]
+  for i=1:length(Ops), j=1:length(d)
+    evals[i][j] = (v[:, j]'*Ops[i].data*v[:, j])[1]
+    if !isapprox(Ops[i].data*v[:, j], evals[i][j]*v[:, j]; atol=atol, rtol=rtol)
+      error("Simultaneous diagonalization failed!")
+    end
+  end
+
+  index = sortperm(real(evals[1][:]))
+  evals_sorted = [real(evals[i][index]) for i=1:length(Ops)]
+  evals_sorted, v[:, index]
+end
+
 end # module

test/test_spectralanalysis.jl

@@ -70,3 +70,24 @@ b = eigenstates_hermitian(H)[1:(Ncutoff+1)]
 for i=1:length(b)
     @test 1-abs(dagger(b[i])*basisstates[i])<1e-12
 end
+
+twospinbasis = spinbasis ⊗ spinbasis
+Sx = full(sum([embed(twospinbasis, i, sx) for i=1:2]))/2.
+Sy = full(sum([embed(twospinbasis, i, sy) for i=1:2]))/2.
+Sz = full(sum([embed(twospinbasis, i, sz) for i=1:2]))/2.
+Ssq = Sx^2 + Sy^2 + Sz^2
+d, v = simdiag([Sz, Ssq])
+@test d[1] == [-1.0, 0, 0, 1.0]
+@test isapprox(d[2],  [2, 0.0, 2, 2])
+@test_throws ErrorException simdiag([Sx, Sy])
+
+threespinbasis = spinbasis ⊗ spinbasis ⊗ spinbasis
+Sx3 = full(sum([embed(threespinbasis, i, sx) for i=1:3])/2.)
+Sy3 = full(sum([embed(threespinbasis, i, sy) for i=1:3])/2.)
+Sz3 = full(sum([embed(threespinbasis, i, sz) for i=1:3])/2.)
+Ssq3 = Sx3^2 + Sy3^2 + Sz3^2
+d3, v3 = simdiag([Ssq3, Sz3])
+dsq3_std = eigvals(full(Ssq3).data)
+@test isapprox(diagm(dsq3_std), v3'*Ssq3.data*v3)
+@test_throws ErrorException simdiag([Sy3, Sz3])
+@test_throws ErrorException simdiag([full(destroy(fockbasis)), full(create(fockbasis))])