Building projection operators efficiently

Question

The projection operator $P$ is defined in quantum mechanics as $$ P=\sum_i^n \vert\psi_i\rangle \langle \psi_i\vert, $$

where $\vert\psi_i\rangle$ is a column vector and $(\langle \psi_i\vert=(\vert\psi_i\rangle)^*)^T$ is a row vector. The $\vert\psi_i\rangle$ that I am using are eigenvectors of a Hermitian matrix.

So here is what I am using:

{energy, states} = Eigensystem[H];

Reenergy = N[Re[energy]];

NormalizedOrderedStates = states[[Ordering[Reenergy]]];

psiket = Map[ConjugateTranspose[{#}] &,NormalizedOrderedStates[[1 ;;n]]]; (*This step turns the states into nx1 matrices*)

psibra = Map[ {#} &, NormalizedOrderedStates[[1 ;;n]]]; (*This step turns the states into 1xn matrices*)

Proj = Sum[psiket[[i]].psibra[[i]], {i, 1, n}];

The problem that I am having is that $H$ can be somewhere from 5000 by 5000 to 50000 by 50000 matrix. The above code is too slow and uses too much memory and in general doesn't feel right. Is there a better way to achieve this?

Answer 1

If I understand it right:

size = 500; H = RandomComplex[1 + I, {size, size}];
n = size; (* set n = size as an example. *)

{energysorted, statessorted} = Eigensystem[H]\[Transpose] // Sort // Transpose;

(bra = statessorted[[;; n]];
 ket = Conjugate@bra;
 myProj = ket\[Transpose] . bra;) // AbsoluteTiming
(* {0.0158231, Null} *)

Key point to understand the code:

1. Check carefully about how Sort works by reading the Details section of its document.

2. Forget about column and row when using Mathematica to deal with matrix.

For comparison:

{energy, states} = Eigensystem[H];

(Reenergy = N[Re[energy]];
  
  NormalizedOrderedStates = states[[Ordering[Reenergy]]];
  
  psiket = Map[ConjugateTranspose[{#}] &, NormalizedOrderedStates[[1 ;; n]]]; 
  
  psibra = Map[{#} &, NormalizedOrderedStates[[1 ;; n]]]; 
  
  Proj = Sum[psiket[[i]] . psibra[[i]], {i, 1, n}];) // AbsoluteTiming
(* {3.67195, Null} *)

As we can see, it's a 200x speed-up.

Because of numeric error, the matrix won't be exactly the same, but:

myProj - Proj // Chop // Flatten // Union
(* {0} *)

BTW, being poor in math, I'm not sure about the reason, but the sorting step seems to be redundant for n = size:

{energy, states} = Eigensystem[H];

(bra = states[[;; n]];
  ket = Conjugate@bra;
  myProj2 = ket\[Transpose] . bra;) // AbsoluteTiming

myProj - myProj2 // Chop // Flatten // Union
(* {0} *)