Update decomposition handling (#364)

* Start purging LAPACK

* Fix typo

* Simplest purge first

* Rework decompositions

* Fix docs build, be consistent in symbol capitalization

* Rerun modified examples

* Fix cotangent sensitivity warning spam

* Split of phase computation in gauge fixing

* Dense defaults

* Gauge fixing updates in tests

* Try to get out of weird sylvester error

* Remove `condition_number`

* Set concrete truncation strategy directly into decomposition alg

* Don't cut between degeneracies

* Avoid Zygote issues `@set` of nested fields

* Actually don't truncate through degeneracies by matching degeneracy with truncation space dimension

* Import `diagview`

* Alternative approach, directly work with `TruncatedAlgorithm`

* Actually set truncation in decomposition algorithm

* Missing import

* Fix trunc setting in SVD wrapper test

* Also don't cut between degeneracies in the eigh wrapper test

* One more missing import

Author: leburgel

Project.toml

@@ -38,7 +38,7 @@ LinearAlgebra = "1"
 LoggingExtras = "1"
 MPSKit = "0.13.9"
 MPSKitModels = "0.4"
-MatrixAlgebraKit = "0.6"
+MatrixAlgebraKit = "0.6.5"
 OhMyThreads = "0.7, 0.8"
 OptimKit = "0.4"
 Printf = "1"

docs/make.jl

@@ -26,6 +26,7 @@ links = InterLinks(
     "MPSKitModels" => "https://quantumkithub.github.io/MPSKitModels.jl/dev/",
     # "Zygote" => "https://fluxml.ai/Zygote.jl/stable/",
     "ChainRulesCore" => "https://juliadiff.org/ChainRulesCore.jl/stable/",
+    "MatrixAlgebraKit" => "https://quantumkithub.github.io/MatrixAlgebraKit.jl/stable/",
 )
 
 # explicitly set math engine

docs/src/examples/3d_ising_partition_function/index.md

@@ -1,16 +1,17 @@
 {
  "cells": [
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "using Markdown #hide"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "# The 3D classical Ising model\n",
     "\n",
@@ -36,25 +37,25 @@
     "spirit as the boundary MPS methods demonstrated in another example.\n",
     "\n",
     "Let's start by making the example deterministic and importing the required packages:"
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "using Random\n",
     "using LinearAlgebra\n",
     "using PEPSKit, TensorKit\n",
     "using KrylovKit, OptimKit, Zygote\n",
     "\n",
     "Random.seed!(81812781144);"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "## Defining the partition function\n",
     "\n",
@@ -65,12 +66,13 @@
     "constituent rank-6 `PEPSKit.PEPOTensor` `O` located at each site of the cubic\n",
     "lattice. To verify our example we will check the magnetization and energy, so we also define\n",
     "the corresponding rank-6 tensors `M` and `E` while we're at it."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "function three_dimensional_classical_ising(; beta, J = 1.0)\n",
     "    K = beta * J\n",
@@ -108,31 +110,30 @@
     "\n",
     "    return TensorMap(o, TMS), TensorMap(m, TMS), TensorMap(e, TMS)\n",
     "end;"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "Let's initialize these tensors at inverse temperature $\\beta=0.2391$, which corresponds to\n",
     "a slightly lower temperature than the critical value $\\beta_c=0.2216544…$"
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "beta = 0.2391\n",
     "O, M, E = three_dimensional_classical_ising(; beta)\n",
     "O isa PEPSKit.PEPOTensor"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "## Contracting the partition function\n",
     "\n",
@@ -190,16 +191,18 @@
     "contraction algorithm we can use to compute the values of these two networks. In addition,\n",
     "we'll specify the specific reverse rule algorithm that will be used to compute the gradient\n",
     "of this cost function."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "boundary_alg = SimultaneousCTMRG(; maxiter = 150, tol = 1.0e-8, verbosity = 1)\n",
     "rrule_alg = EigSolver(;\n",
-    "    solver_alg = KrylovKit.Arnoldi(; maxiter = 30, tol = 1.0e-6, eager = true), iterscheme = :diffgauge\n",
+    "    solver_alg = KrylovKit.Arnoldi(; maxiter = 30, tol = 1.0e-6, eager = true),\n",
+    "    iterscheme = :fixed,\n",
     ")\n",
     "T = InfinitePEPO(O)\n",
     "\n",
@@ -240,12 +243,11 @@
     "    g = only(gs)\n",
     "    return E, g\n",
     "end;"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "There are a few things to note about this cost function definition. Since we will pass it to\n",
     "the `OptimKit.optimize`, we require it to return both our cost function and the\n",
@@ -286,12 +288,13 @@
     "them where the only difference is that we have to pass along an extra environment since our\n",
     "cost function requires two distinct contractions as opposed to the setting of Hamiltonian\n",
     "PEPS optimization which only requires a double-layer contraction."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "function pepo_retract((peps, env_double_layer, env_triple_layer), η, α)\n",
     "    (peps´, env_double_layer´), ξ = PEPSKit.peps_retract((peps, env_double_layer), η, α)\n",
@@ -309,24 +312,24 @@
     "        ξ, (peps, env_double_layer), η, α, (peps´, env_double_layer´)\n",
     "    )\n",
     "end;"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "### Finding the fixed point\n",
     "\n",
     "All that is left then is to specify the virtual spaces of the PEPS and the two environments,\n",
     "initialize them in the appropriate way, choose an optimization algortithm and call the\n",
     "`optimize` function from OptimKit.jl to get our desired PEPS fixed point."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "Vpeps = ℂ^2\n",
     "Venv = ℂ^12\n",
@@ -345,41 +348,41 @@
     "    retract = pepo_retract,\n",
     "    (transport!) = (pepo_transport!),\n",
     ");"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "### Verifying the result\n",
     "\n",
     "Having found the fixed point, we have essentially contracted the entire partition function\n",
     "and we can start computing observables. The free energy per site for example is just given by\n",
     "the final value of the cost function we have just optimized."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "@show f"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "As another check, we can compute the magnetization per site and compare it to a [reference\n",
     "value obtaind through Monte-Carlo simulations](@cite hasenbusch_monte_2001)."
-   ],
-   "metadata": {}
+   ]
   },
   {
-   "outputs": [],
    "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "n3_final = InfiniteSquareNetwork(psi_final, T)\n",
     "num = PEPSKit.contract_local_tensor((1, 1, 1), M, n3_final, env3_final)\n",
@@ -389,33 +392,31 @@
     "m_ref = 0.667162\n",
     "\n",
     "@show abs(m - m_ref)"
-   ],
-   "metadata": {},
-   "execution_count": null
+   ]
   },
   {
    "cell_type": "markdown",
+   "metadata": {},
    "source": [
     "---\n",
     "\n",
     "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
-   ],
-   "metadata": {}
+   ]
   }
  ],
- "nbformat_minor": 3,
  "metadata": {
+  "kernelspec": {
+   "display_name": "Julia 1.12.5",
+   "language": "julia",
+   "name": "julia-1.12"
+  },
   "language_info": {
    "file_extension": ".jl",
    "mimetype": "application/julia",
    "name": "julia",
    "version": "1.12.5"
-  },
-  "kernelspec": {
-   "name": "julia-1.12",
-   "display_name": "Julia 1.12.5",
-   "language": "julia"
   }
  },
- "nbformat": 4
+ "nbformat": 4,
+ "nbformat_minor": 3
 }