-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
19 changed files
with
592 additions
and
176 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,234 @@ | ||
| { | ||
| "cells": [ | ||
| { | ||
| "cell_type": "raw", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "---\n", | ||
| "title: \"`primate` usage - quickstart\"\n", | ||
| "---" | ||
| ], | ||
| "id": "032be247" | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "Below is a quick introduction to `primate`. For more introductory material, theor\n" | ||
| ], | ||
| "id": "d2a9b9ea" | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "#| echo: false\n", | ||
| "#| output: false\n", | ||
| "from bokeh.plotting import figure, show\n", | ||
| "from bokeh.io import output_notebook\n", | ||
| "output_notebook()\n", | ||
| "import numpy as np\n", | ||
| "np.random.seed(1234)" | ||
| ], | ||
| "id": "d88964cc", | ||
| "execution_count": null, | ||
| "outputs": [] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "To do trace estimation, use functions in the `trace` module: " | ||
| ], | ||
| "id": "9bb3cfa6" | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "import primate.trace as TR\n", | ||
| "from primate.random import symmetric\n", | ||
| "A = symmetric(150) ## random positive-definite matrix \n", | ||
| "\n", | ||
| "print(f\"Actual trace: {A.trace():6f}\") ## Actual trace\n", | ||
| "print(f\"Girard-Hutch: {TR.hutch(A):6f}\") ## Monte-carlo tyoe estimator\n", | ||
| "print(f\"XTrace: {TR.xtrace(A):6f}\") ## Epperly's algorithm" | ||
| ], | ||
| "id": "ba7ff22b", | ||
| "execution_count": null, | ||
| "outputs": [] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "For matrix functions, you can either construct a `LinearOperator` directly via the `matrix_function` API, or supply a string to the parameter `fun` describing the spectral function to apply. For example, one might compute the log-determinant as follows:\n" | ||
| ], | ||
| "id": "2fc4b61b" | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "from primate.operator import matrix_function\n", | ||
| "M = matrix_function(A, fun=\"log\")\n", | ||
| "\n", | ||
| "ew = np.linalg.eigvalsh(A)\n", | ||
| "print(f\"logdet(A): {np.sum(np.log(ew)):6f}\")\n", | ||
| "print(f\"GR approx: {TR.hutch(M):6f}\")\n", | ||
| "print(f\"XTrace: {TR.xtrace(M):6f}\")\n", | ||
| "\n", | ||
| "## Equivalently could've used: \n", | ||
| "## M = matrix_function(A, fun=np.log)" | ||
| ], | ||
| "id": "e07ecdd1", | ||
| "execution_count": null, | ||
| "outputs": [] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "Note in the above example you can supply to `fun` either string describing a built-in spectral function or an arbitrary `Callable`. The former is preferred when possible, as function evaluations will generally be faster and `hutch` can also be parallelized. Multi-threaded execution of e.g. `hutch` with arbitrary functions is not currently allowed due to the GIL, though there are options available, see [the integration docs](../advanced/cpp_integration.qmd) for more details. \n", | ||
| "\n", | ||
| "For 'plain' operators, `XTrace` should recover the exact trace (up to roundoff error). For matrix functions $f(A)$, there will be some inherent inaccuracy as the underlying matrix-vector multiplication is approximated with the Lanczos method. \n", | ||
| "\n", | ||
| "\n", | ||
| "\n", | ||
| "\n", | ||
| "\n", | ||
| "\n", | ||
| "\n", | ||
| "\n", | ||
| "In general, the amount of accuracy depends both on the Lanczos parameters and the type of matrix function. Spectral functions that are difficult or impossible to approximate via low-degree polynomials, for example, may suffer more from inaccuracy issues than otherwise. For example, consider the example below that computes that rank:\n" | ||
| ], | ||
| "id": "8fe4ef85" | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "## Make a rank-deficient operator\n", | ||
| "ew = np.sort(ew)\n", | ||
| "ew[:30] = 0.0\n", | ||
| "A = symmetric(150, ew = ew, pd = False)\n", | ||
| "M = matrix_function(A, fun=np.sign)\n", | ||
| "\n", | ||
| "print(f\"Rank: {np.linalg.matrix_rank(A)}\")\n", | ||
| "print(f\"GR approx: {TR.hutch(M)}\")\n", | ||
| "print(f\"XTrace: {TR.xtrace(M)}\")" | ||
| ], | ||
| "id": "054e65ba", | ||
| "execution_count": null, | ||
| "outputs": [] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "This is not so much a fault of `hutch` or `xtrace` as much as it is the choice of approximation and Lanczos parameters. The `sign` function has a discontinuity at 0, is not smooth, and is difficult to approximate with low-degree polynomials. \n", | ||
| "One workaround to handle this issue is relax the sign function with a low-degree \"soft-sign\" function: \n", | ||
| "$$ \\mathcal{S}_\\lambda(x) = \\sum\\limits_{i=0}^q \\left( x(1 - x^2)^i \\prod_{j=1}^i \\frac{2j - 1}{2j} \\right)$$\n", | ||
| "\n", | ||
| "Visually, the soft-sign function looks like this: " | ||
| ], | ||
| "id": "931412fd" | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "from primate.special import soft_sign, figure_fun\n", | ||
| "show(figure_fun(\"smoothstep\"))" | ||
| ], | ||
| "id": "8a3f3498", | ||
| "execution_count": null, | ||
| "outputs": [] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "It's been shown that there is a low degree polynomial $p^\\ast$ that uniformly approximates $\\mathcal{S}_\\lambda$ up to a small error on the interval $[-1,1]$. Since `matrix_function` uses a low degree Krylov subspace to approximate the action $v \\mapsto f(A)v$, one way to improve the accuracy of rank estimation is to replace $\\mathrm{sign} \\mapsto \\mathcal{S}_{\\lambda}$ for some choice of $q \\in \\mathbb{Z}_+$ (this function is available in `primate` under the name `soft_sign`):\n" | ||
| ], | ||
| "id": "4c10e2d4" | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "from primate.special import soft_sign\n", | ||
| "for q in range(0, 50, 5):\n", | ||
| " M = matrix_function(A, fun=soft_sign(q=q))\n", | ||
| " print(f\"XTrace S(A) for q={q}: {TR.xtrace(M):6f}\")" | ||
| ], | ||
| "id": "b50e55c9", | ||
| "execution_count": null, | ||
| "outputs": [] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "If the type of operator `A` is known to typically have a large [spectral gap](), another option is to compute the _numerical rank_ by thresholding values above some fixed value $\\lambda_{\\text{min}}$. This is equivalent to applying the following spectral function:\n", | ||
| "\n", | ||
| "$$ S_{\\lambda_{\\text{min}}}(x) = \n", | ||
| "\\begin{cases} \n", | ||
| "1 & \\text{ if } x \\geq \\lambda_{\\text{min}} \\\\\n", | ||
| "0 & \\text{ otherwise }\n", | ||
| "\\end{cases}\n", | ||
| "$$\n", | ||
| "\n", | ||
| "In the above example, the optimal cutoff $\\lambda_{\\text{min}}$ is given by the smallest non-zero eigenvalue. Since the trace estimators all stochastic to some degree, we set the cutoff to slightly less than this value: \n" | ||
| ], | ||
| "id": "07810e39" | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "lambda_min = min(ew[ew != 0.0])\n", | ||
| "print(f\"Smallest non-zero eigenvalue: {lambda_min:.6f}\")\n", | ||
| "\n", | ||
| "step = lambda x: 1 if x > (lambda_min * 0.90) else 0\n", | ||
| "M = matrix_function(A, fun=step, deg=50)\n", | ||
| "print(f\"XTrace S_t(A) for t={lambda_min*0.90:.4f}: {TR.xtrace(M):6f}\")" | ||
| ], | ||
| "id": "ddb4dea7", | ||
| "execution_count": null, | ||
| "outputs": [] | ||
| }, | ||
| { | ||
| "cell_type": "markdown", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "Indeed, this works! Of course, here we've used the fact that we know the optimal cutoff value, but this can also be estimated with the `lanczos` method itself. \n" | ||
| ], | ||
| "id": "e82a40b8" | ||
| }, | ||
| { | ||
| "cell_type": "code", | ||
| "metadata": {}, | ||
| "source": [ | ||
| "from primate.diagonalize import lanczos\n", | ||
| "from scipy.linalg import eigvalsh_tridiagonal\n", | ||
| "a,b = lanczos(A)\n", | ||
| "rr = eigvalsh_tridiagonal(a,b) # Rayleigh-Ritz values\n", | ||
| "tol = 10 * np.finfo(A.dtype).resolution\n", | ||
| "print(f\"Approx. cutoff: {np.min(rr[rr > tol]):.6f}\")" | ||
| ], | ||
| "id": "b1ba0b78", | ||
| "execution_count": null, | ||
| "outputs": [] | ||
| } | ||
| ], | ||
| "metadata": { | ||
| "kernelspec": { | ||
| "name": "python3", | ||
| "language": "python", | ||
| "display_name": "Python 3 (ipykernel)" | ||
| } | ||
| }, | ||
| "nbformat": 4, | ||
| "nbformat_minor": 5 | ||
| } |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.