jcmgray
.H property on MPO performs only complex conjugation, not full Hermitian conjugate (dagger)
What is your issue?
I am using Quimb with the PyTorch backend to work with Matrix Product Operators (MPOs). I noticed that the .H property on an MPO does not perform the full Hermitian conjugate (dagger) operation as expected. Instead, it only applies the complex conjugate to each local tensor, without transposing the physical indices or reversing the order of the MPO tensors.
Environment:
quimb version: 1.11.0
Python version: 3.10
Operating system: Linux
Installation method: conda environment
import quimb as qu
import quimb.tensor as qtn
import torch
def print_tns(mpo):
# Extract the raw arrays and a skeleton of the TN
params, skeleton = qtn.pack(mpo)
for p in params:
print(params[p])
def main():
L = 2
mpo_bond_dim = 1
mpo = qtn.MatrixProductOperator.from_fill_fn(
fill_fn=lambda shape: torch.rand(*shape, dtype=torch.complex64),
L=L,
bond_dim=mpo_bond_dim,
phys_dim=4,
tags='MyMPO',
upper_ind_id='UP{}',
lower_ind_id='LO{}'
)
print(mpo)
print_tns(mpo)
print('mpo.H')
print_tns(mpo.H)
if __name__ == "__main__":
main()
Observed output:
MatrixProductOperator([
Tensor(shape=(1, 4, 4), inds=('_4ae00cAAAAB', 'UP0', 'LO0'), tags=oset(['MyMPO', 'I0'])),
Tensor(shape=(1, 4, 4), inds=('_4ae00cAAAAB', 'UP1', 'LO1'), tags=oset(['MyMPO', 'I1'])),
], tensors=2, indices=5, L=2, max_bond=1)
tensor([[[0.7632+0.7197j, 0.4769+0.2170j, 0.5821+0.5419j, 0.2200+0.3946j],
[0.5214+0.4283j, 0.4174+0.7914j, 0.7996+0.0899j, 0.5664+0.7025j],
[0.8373+0.9964j, 0.2938+0.6747j, 0.7179+0.4150j, 0.4569+0.5680j],
[0.0818+0.4894j, 0.6299+0.1808j, 0.0332+0.7765j, 0.5294+0.2070j]]])
tensor([[[0.7724+0.9414j, 0.8078+0.5182j, 0.4615+0.6405j, 0.0179+0.4112j],
[0.0736+0.3292j, 0.5040+0.0107j, 0.2087+0.4135j, 0.5105+0.5459j],
[0.7177+0.2983j, 0.4167+0.7490j, 0.0522+0.1469j, 0.0736+0.3903j],
[0.2155+0.9970j, 0.6950+0.6002j, 0.6959+0.1622j, 0.4879+0.9362j]]])
mpo.H
tensor([[[0.7632-0.7197j, 0.4769-0.2170j, 0.5821-0.5419j, 0.2200-0.3946j],
[0.5214-0.4283j, 0.4174-0.7914j, 0.7996-0.0899j, 0.5664-0.7025j],
[0.8373-0.9964j, 0.2938-0.6747j, 0.7179-0.4150j, 0.4569-0.5680j],
[0.0818-0.4894j, 0.6299-0.1808j, 0.0332-0.7765j, 0.5294-0.2070j]]])
tensor([[[0.7724-0.9414j, 0.8078-0.5182j, 0.4615-0.6405j, 0.0179-0.4112j],
[0.0736-0.3292j, 0.5040-0.0107j, 0.2087-0.4135j, 0.5105-0.5459j],
[0.7177-0.2983j, 0.4167-0.7490j, 0.0522-0.1469j, 0.0736-0.3903j],
[0.2155-0.9970j, 0.6950-0.6002j, 0.6959-0.1622j, 0.4879-0.9362j]]])
Clearly, .H only applies complex conjugation and does not transpose the tensors or reverse their order.
This behavior seems like a bug or at least a misleading implementation since the .H property name suggests the full Hermitian conjugate.
The Quimb + PyTorch tutorial examples appear to imply .H behaves as a full dagger, which adds to the confusion.
It would be helpful if .H either performed the full dagger (conjugate transpose + reverse order) or if the documentation clarified its behavior explicitly.
Also, guidance on how to correctly implement the full dagger operation for MPOs in Quimb with the PyTorch backend would be appreciated.
If useful, I can help prepare a minimal reproducible example or propose a fix/PR for this behavior.
This is a bit confusing but intentional, the way indices works means that whether something has been 'transposed' in the matrix sense depends on how you compose it with other objects, for example mpo.H @ mpo indeed acts like it has been transposed. One just has to learn a bit how TNs are composed in quimb, they don't act like matrices.
To form the dagger for calling apply, one takes the conjugate, and then just needs to swap the lower and upper ind labels - one doesn't need to do anything to the arrays (including reversing any order, that would be like permuting the physical dimensions). Here's an implementation that will work for any operator like TN:
import quimb.tensor as qtn
def dagger(op):
ixmap = {
op.lower_ind(site): op.upper_ind(site)
for site in op.sites
} | {
op.upper_ind(site): op.lower_ind(site)
for site in op.sites
}
return op.conj().reindex_(ixmap)
Such a method could be added TensorNetworkGenOperator, which MPOs and PEPOs inherit from, if it was useful, but it would need to well documented that a.dagger() @ a != a.H @ a and does not work like a inner product.
