4. Linear Algebra

4.1. Dense / full eigen decomposition

Currently full decompositions use numpy. They are as follows:

• eig()

• eigh()

• eigvals()

• eigvalsh()

• eigvecs()

• eigvecsh()

• eigensystem()

• svd()

4.1.1. Automatic Exploitation of Symmetries via Blocking

Sometimes symmetries present themselves in an operator as, after a certain permutation, diagonal blocks. This fact can be exploited to speed up full eigen-decompositions as each block can be decomposed individually. For example, imagine we want to supply the eigen-decomposition of the Heisenberg Hamiltonian to an Evolution object to perform long-time, exact, dynamics:

import quimb as qu

H = qu.ham_heis(12)

%%timeit
el, ev = qu.eigh(H)

9.17 s ± 850 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

This is already quite slow, and will not be sustainable for larger lengths. The Heisenberg Hamiltonian however has an abelian \(Z_2\) subgroup symmetry, conserved Z-spin, that manifests itself in the computational basis. Thus if instead we specify autoblock=True:

%%timeit
el, ev = qu.eigh(H, autoblock=True)

466 ms ± 63.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

We can see a drastic speedup (helped by the fact that the whole algorithm uses numba). To work more fundamentally with the indentified blocks see the compute_blocks() function.

4.2. Partial eigen decomposition

Partial decompositions are mostly just specified by supplying the k kwarg to the above functions. These also take a backend argument which can be one of:

• 'scipy': is generally reliable

• 'numpy': can be faster for small or dense problems

• 'lobpcg': useful for fast, low accruacy generalized eigenproblems (like periodic DMRG)

• 'slepc': Usually the fastest for large problems, with many options. Will either spawn MPI workers or should be used in syncro mode.

• 'slepc-nompi': like 'slepc', but performs computation in the single, main process.

• 'AUTO' - choose a good backend, the default.

The possible partical decompositions are:

• eig() with (k > 0)

• eigh() with (k > 0)

• eigvals() with (k > 0)

• eigvalsh() with (k > 0)

• eigenvectors() with (k > 0)

• eigvecs() with (k > 0)

• eigvecsh() with (k > 0)

• groundstate()

• groundenergy()

• eigensystem_partial()

• svds()

So for example the groundstate() function for a Hamiltonian H is an alias to:

psi = qu.eigvecsh(H, k=1, which='SA')

[find eigenvectors, Hermitian operator (h post-fix), get k=1 eigenstate, and target the ‘(s)mallest (a)lgebraic’ eigenvalue].

4.3. Interior eigen-solving

SLEPc is highly recommended for performing these using ‘shift-invert’. See the following functions:

• eigh(..., k=k, sigma=x) with k > 0 etc., or

• eigh_window()

• eigvalsh_window()

• eigvecsh_window()

With the last three allowing the specification of a window relative to the total spectrum of the operator.

4.4. Fast Randomized Linear Algebra

quimb has an implementation of a fast randomized SVD - rsvd() - that can be significantly quicker than svd() or svds(), especially for large k. This might be useful for e.g. tensor network linear operator decompositions. It can perform the SVD rank-adaptively, which allows the efficient estimation of an operator’s rank, see estimate_rank().