Building Projection Operators Onto Subspaces

Question

Intro

In quantum mechanics one uses projection operators in many problems / derivations.

At the problem at hand I need to build projection operators (matrices) onto a subspace of my Hilbert space. I have a 16 dimensional Hilbert space and I wish to perform a spectral decomposition (in the Hamiltonian's spectrum) on arbitrary operators using projections in the following manner. Let us assume that I have an operator $A$, then I could define:

$$ A(\omega)=\sum_{\epsilon'-\epsilon=\omega} P(\epsilon) A P(\epsilon') $$

Then $A(\omega)$ is a projection of $A$,and we have the completeness relation $\sum_\omega A(\omega)=A$.

For clarity the sum runs over all energies of the system (eigenvalues of the Hamiltonian), such that the difference between the two is equal to $\omega$. the projection $P(\epsilon)$ stands for the projection over the subspace of Hilbert space which belongs to energy $\epsilon$.

For those interested my motivation to do this decomposition follows from chap.3 of The Theory of Open Quantum Systems. Another example may be found free online here.

Given that I'm able to diagonalize my Hamiltonian (or else this discussion is futile), I can trivially rotate $A$ to the eigenbasis of of the Hamiltonian in which the projection operators reduce to diagonal matrices with 0's and 1's on the diagonal which makes it easier to build them.

My Questions:

1. Building $P(\epsilon)$

As explained, naturally I must diagonalize the Hamiltonian first. However I'm unable to do that symbolically, but only numerically, how would I build the operators $P(\epsilon)$? To emphasize I might need to do this many times (iteratively), so it cannot be done by hand. Per iteration I need:

1. To diagonalize numerically (trivial)
2. Infer which which energies do I have. As my spectrum is necessarily degenerate this may not be so trivial. Since I diagonalize numerically how do I decide if two energies are close enough to be classified as degenerate?
3. Build the operators for each level. Can it be done using fast built-in function, or do I have to iterate using For loop which will be very inefficient?

2. Summation

In order to do the summation I come up with the following scheme:

1. Find out which energy differences are possible in the system which isn't too hard as I have the energies.
2. Preform a double summation with a KroneckerDelta.

First,I'm not sure that a double summation is the most efficient way, is there a better one? Secondary, how will the KroneckerDelta preform with numerical values? Will KroneckerDelta [Chop[e1-e2-omega]], where e1 and e2 are my supposed summation indices work?

My Attempt

Here's my poor man solution for the problem, I love to get responses on how to improve it / save some computation time or steps. This solution can be incorporated into a function of course, but I don't do it here for the sake of readability.

Given a Hamiltonian (Hermitian matrix) h according to which eigenvalue subspaces one wants to decompose a second operator (matrix) Op,written in the same basis, one first need to decompose h.

{e, v} = Chop[Eigensystem[h]];
uniqueE = DeleteDuplicates[e];

Using the list of unique energies (eigenvalues) on may build the projection operators. These are most easily built in the eigenbasis of the Hamiltonian as in this basis these are just diagonal matrices, with zeros and ones one there diagonal:

projectionOP = ParallelMap[{#, 
  SparseArray[{i_, i_} :> Boole[e[[i]] == #], Dimensions[h]]} &,uniqueE];

This returns a list of pairs, where each pair consists of one of the uniqueE (eigenvalue of h), and a diagonal matrices with entries along its diagonal of 0 and 1, and we have 1 where e[[i]] is equal to the relevant element in uniqueE.

Since I built the matrices in projectionOP in the eigenbasis of h I'll rotate Op to the same basis (though it could be done the there way around):

rotatedOp = v.Op.Transpose[v]

Now I get to the part of my implementation with which I'm not so happy, as I feel it could be much simpler. First I find all of the projections, i.e. all of the possible multiplications which I wrote on top regardless of the value of $\omega$ written there:

AllProjections = 
Flatten[Outer[{#2[[1]] - #1[[1]], #1[[2]]. rotatedOp.#2[[2]]} &, 
 projectionOP, projectionOP, 1], 1];

This returns a list of pairs where each pair consists or a certain $\omega = \epsilon'-\epsilon$, or a in the MMA syntax uniqueE[[i]]-uniqueE[[j]] for some i and j.

The next step may not be a necessary one but it is convenient. Some of the matrix multiplications may vanish depending on their specific form. For instance in my case the matrix Op is rather sparse, and since the projection operators are also sparse, most of the matrices in the above expression are all zeros. ThusI remove them:

nonZeroProjections = 
Select[AllProjections, #[[2]] != ConstantArray[0, Dimensions[h]] &];

Lastly I wish to preform the summation (as there might be several contributions per $\omega$). Therefore I GatherBy the different pairs according to their value of $\omega$, and sum:

DecomposedOp = 
Map[{#[[1, 1]], Plus @@ #[[All, 2]]} &, GatherBy[nonZeroProjections, First]];

The completeness relation (I think the set of projection operators forms a complete basis) may be verified easily by:

(Plus @@ DecomposedOp[[All, 2]]) - rotatedOp == ConstantArray[0, Dimensions[h]]

Remaining Issues

1. The matrices projectionOp have an amount of 1's on there diagonal which is equal to the degeneracy of each eigenvalue. This part of the implementation is prone to error if due to numeric issues two degenerate eigenvalues don't seem degenerate. I checked my code with a degenerate case and it seems to work there (I know how many unique eigenvalues to expect there). How do I ensure correct behavior of this part?

2. As can be seen after the matrices projectionOp are found, the summation still takes several steps to preform. Any better way to do this?

Answer 1

I presume that you use the Euclidean scalarproduct for diagonalizing the Hamiltonian. Otherwise you would use the generalized eigensystem facilities of Eigensystem or a CholeskyDecomposition of the (inverse) of the Gram matrix.

Let's generate some example data.

H1 = RandomReal[{-1, 1}, {160, 160}];
H1 = Transpose[H1].H1;
H = ArrayFlatten[{
    {H1, 0., 0., 0.},
    {0., H1, 0., 0.},
    {0., 0., H1, 0.},
    {0., 0., 0., H1 + 0.000000001}
    }];
A = RandomReal[{-1, 1}, Dimensions[H]];

The interesting parts starts here. I use ClusteringComponents to find clusters within the eigenvalues and their differences. This should make it a bit more robust.

{lambda, U} = Eigensystem[H];
eigclusters = GroupBy[
   Transpose[{ClusteringComponents[lambda], Range[Length[H]]}],
   First -> Last
   ];
P = Association[
   Map[x \[Function] Mean[lambda[[x]]] -> Transpose[U[[x]]].U[[x]], 
    Values[eigclusters]]];

diffs = Flatten[Outer[Plus, Keys[P], -Keys[P]], 1];
pos = Flatten[Outer[List, Range[Length[P]], Range[Length[P]]], 1];
diffclusters = GroupBy[
   Transpose[{ClusteringComponents[diffs], Range[Length[diffs]]}],
   First -> Last
   ];
prA = Association[
   Map[
    x \[Function] Rule[
      Mean[diffs[[x]]],
      Total[Map[y \[Function] P[[y[[1]]]].A.P[[y[[2]]]], pos[[x]]]
       ]
      ],
    Values[diffclusters]
    ]
   ];

prA is now an association that maps energy differences to the according projections of A. You may check with Max[Abs[Total[prA] - A]] that the decomposition is complete.

Note however, that the eigenspaces are instable under perturbation. This is a mathematical issue, not any that is related to Mathematica.

Moreover, ClusteringComponents really slows things down. But there are plenty of opportunities to configure ClusteringComponents. Surely, you can tune it to your needs. See the documentation on ClusteringComponents.

Edit: A slight speedup when dealing with larger matrices can be obtained by leaving the projectors in factorized form:

{lambda, U} = Eigensystem[H];
eigclusters = 
  GroupBy[Transpose[{ClusteringComponents[lambda], Range[Length[H]]}],
    First -> Last];
V = Association[
   Map[x \[Function] Mean[lambda[[x]]] -> U[[x]], 
    Values[eigclusters]]];

diffs = Flatten[Outer[Plus, Keys[P], -Keys[P]], 1];
pos = Flatten[Outer[List, Range[Length[P]], Range[Length[P]]], 1];
diffclusters = 
  GroupBy[Transpose[{ClusteringComponents[diffs], 
     Range[Length[diffs]]}], First -> Last];
prA2 = Association[Map[x \[Function] Rule[
      Mean[diffs[[x]]],
      Total[
       Map[
        y \[Function] With[{U1 = V[[y[[1]]]], U2 = V[[y[[2]]]]},
          Transpose[U1].(U1.A.Transpose[U2]).U2
          ],
        pos[[x]]
        ]
       ]
      ],
    Values[diffclusters]
    ]
   ];

But the cost of ClusteringComponents is still dominant.