Cirrus changes

.cirrus.yml

@@ -4,18 +4,17 @@ env:
 linux_task:
   container:
     image: quay.io/pypa/manylinux_2_28_x86_64
-  setup_python_script: |
-    export PATH=$PATH:/opt/python/cp310-cp310/bin/
-    alias python=/opt/python/cp310-cp310/bin/python
+  only_if: changesInclude('.cirrus.yml', '**.{h,cpp,py}') 
+  env:
+    CC: clang 
+    CXX: clang++
+    PATH: ${PATH}:/opt/python/cp310-cp310/bin
   pip_cache:
     folder: ~/.cache/pip
     fingerprint_script: echo $PYTHON_VERSION
-  environment:
-    CC: clang 
-    CXX: clang++
   before_script: |
     python -m pip install --upgrade pip
-    python -m pip install pytest pytest-cov pytest-benchmark
+    python -m pip install build pytest pytest-cov pytest-benchmark
   dependencies_script: |
     bash tools/cibw_linux.sh
   build_script: |
@@ -25,7 +24,7 @@ linux_task:
   uninstall_script: |
     python -m pip uninstall primate --yes
   wheel_script: |
-    python -m pip wheel . --no-deps --wheel-dir dist/ --verbose
+    python -m pip build
   install_script: |
     python -m pip install dist/primate*.whl
   test_post_script: |
@@ -35,18 +34,22 @@ linux_task:
 windows_task:
   windows_container:
     image: cirrusci/windowsservercore:2019
+  only_if: changesInclude('.cirrus.yml', '**.{h,cpp,py}') 
   setup_python_script: |
     choco install -y python311
   dependencies_script: |
     bash tools/cibw_windows.sh
+  before_script: |
+    python -m pip install --upgrade pip
+    python -m pip install build pytest pytest-cov pytest-benchmark
   build_script: |
     python -m pip install . --verbose
   test_script: |
     python -m pytest tests/ --cov=primate --benchmark-skip
   uninstall_script: |
     python -m pip uninstall primate -y
   wheel_script: |
-    python -m pip wheel . --no-deps --wheel-dir dist/ --verbose
+    python -m build 
   install_script: |
     python -m pip install dist/primate*.whl
   test_post_script: |

docs/basic/integration.html

@@ -345,14 +345,14 @@ <h1 class="title">Integration</h1>
 </header>
 
 
-<p><code>primate</code> supports a variety of matrix-types of the box, including numpy <code>ndarray</code>’s, compressed <a href="https://docs.scipy.org/doc/scipy/reference/sparse.html">sparse matrices</a> (a lá SciPy), and <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.LinearOperator.html">LinearOperators</a>. Note the last option enables the use of <em>matrix free</em> operators.</p>
+<p><code>primate</code> supports a variety of matrix-types of the box, including numpy <code>ndarray</code>’s, compressed <a href="https://docs.scipy.org/doc/scipy/reference/sparse.html">sparse matrices</a> (a lá SciPy), and <a href="https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.LinearOperator.html">LinearOperators</a>—the latter enables the use of <em>matrix free</em> operators.</p>
 <p>More generally, the basic requirements for any operator <code>A</code> to be used with e.g.&nbsp;the <code>Lanczos</code> method in <code>primate</code> are:</p>
 <ol type="1">
 <li>A method <code>A.matvec(input: ndarray) -&gt; ndarray</code> implementing <span class="math inline">v \mapsto Av</span></li>
 <li>A method <code>A.shape() -&gt; tuple[int, int]</code> giving the output/input dimensions of <span class="math inline">A</span></li>
 </ol>
 <!-- Note that in the matrix setting, `A.shape()` yields `(A.nrow(), A.ncol())`, so existing matrix implementations like [ndarray](https://numpy.org/doc/stable/reference/generated/numpy.ndarray.html)'s and [sparray](https://docs.scipy.org/doc/scipy/reference/sparse.html#sparse-array-classes)'s natively support the interface. -->
-<p>These semantics extend to the C++ side as well via <em>C++20 concepts</em>—see the <a href="advanced/cpp_integration.qmd">C++ integration guide</a>. If you’re using pybind11 and you want to extend <code>primate</code>’s Python API to work natively with linear operator implemented in C++, see the <a href="advanced/pybind11_integration.qmd">pybind11 integration guide</a>.</p>
+<p>These semantics extend to the C++ side as well via <em>C++20 concepts</em>—see the <a href="../advanced/cpp_integration.html">C++ integration guide</a>. If you’re using pybind11 and you want to extend <code>primate</code>’s Python API to work natively with linear operator implemented in C++, see the <a href="../advanced/pybind11_integration.html">pybind11 integration guide</a>.</p>
 
 
 

docs/search.json

@@ -412,7 +412,7 @@
     "href": "basic/integration.html",
     "title": "Integration",
     "section": "",
-    "text": "primate supports a variety of matrix-types of the box, including numpy ndarray’s, compressed sparse matrices (a lá SciPy), and LinearOperators. Note the last option enables the use of matrix free operators.\nMore generally, the basic requirements for any operator A to be used with e.g. the Lanczos method in primate are:\n\nA method A.matvec(input: ndarray) -> ndarray implementing v \\mapsto Av\nA method A.shape() -> tuple[int, int] giving the output/input dimensions of A\n\n\nThese semantics extend to the C++ side as well via C++20 concepts—see the C++ integration guide. If you’re using pybind11 and you want to extend primate’s Python API to work natively with linear operator implemented in C++, see the pybind11 integration guide.",
+    "text": "primate supports a variety of matrix-types of the box, including numpy ndarray’s, compressed sparse matrices (a lá SciPy), and LinearOperators—the latter enables the use of matrix free operators.\nMore generally, the basic requirements for any operator A to be used with e.g. the Lanczos method in primate are:\n\nA method A.matvec(input: ndarray) -> ndarray implementing v \\mapsto Av\nA method A.shape() -> tuple[int, int] giving the output/input dimensions of A\n\n\nThese semantics extend to the C++ side as well via C++20 concepts—see the C++ integration guide. If you’re using pybind11 and you want to extend primate’s Python API to work natively with linear operator implemented in C++, see the pybind11 integration guide.",
     "crumbs": [
       "Basics",
       "Integration"
@@ -423,7 +423,7 @@
     "href": "theory/lanczos.html",
     "title": "The Lanczos method",
     "section": "",
-    "text": "In 1950, Cornelius Lanczos proposed the method of minimized iterations—now known as the Lanczos method—enabling the world to iteratively expose the spectrum of a symmetric A via tridiagonalization. Despite its age, the Lanczos iteration remains the standard1 algorithm both for computing eigensets and for solving linear systems in the large-scale regime.\nIn this post, I’ll describe the basic theory behind the Lanczos method, with a focus on its time and space complexity. In general, the Lanczos method is intrinsically connected to other mathematical constructions, e.g. orthogonal polynomials, gaussian quadrature, conjugate gradient—none of these will be discussed for now.",
+    "text": "In 1950, Cornelius Lanczos proposed the method of minimized iterations—now known as the Lanczos method—enabling the world to iteratively expose the spectra of symmetric operators A \\in \\mathbb{S}_n via tridiagonalization. Despite its age, the Lanczos iteration remains the standard algorithm1 both for computing eigensets and for solving linear systems in the large-scale regime. Listed among the top 10 most influential algorithms of the 20th century in an IEEE guest editorial,\nprimate offers a matrix-free\nIn this post, I’ll describe the basic theory behind the Lanczos method, with a focus on its time and space complexity. In general, the Lanczos method is intrinsically connected to other mathematical constructions, e.g. orthogonal polynomials, gaussian quadrature, conjugate gradient—none of these will be discussed for now.",
     "crumbs": [
       "Theory",
       "Lanczos"
@@ -542,7 +542,7 @@
     "href": "theory/lanczos.html#quadratic-time-and-linear-space-how",
     "title": "The Lanczos method",
     "section": "Quadratic time and linear space? How?",
-    "text": "Quadratic time and linear space? How?\nUnless you know the tricks, its not obvious at all decompositions above takes just O(n^2) time and O(n) space to obtain. So… how does the complexity argument play out?\nIf you squint hard enough, you can deduce that since A Q_j = Q_j T_j + \\beta_{j+1} q_{j+1} e_{j}^T, every symmetric A \\in \\mathbb{R}^{n \\times n} expanded this way admits a three-term recurrence: \n\\begin{align*}\nA q_j &= \\beta_{j\\text{-}1} q_{j\\text{-}1} + \\alpha_j q_j + \\beta_j q_{j+1} \\\\\n\\Leftrightarrow \\beta_{j} \\, q_{j+1} &= A q_j - \\alpha_j \\, q_j - \\beta_{j\\text{-}1} \\, q_{j\\text{-}1}  \n\\end{align*}\n\nThe equation above is a variable-coefficient second-order linear difference equation, and it is known such equations have unique solutions: \n\\alpha_j = q_j^T A q_j, \\;\\; \\beta_j = \\lVert r_j \\rVert_2, \\;\\; q_{j+1} = r_j / \\beta_j\n\n\n\\text{where  } r_j = (A - \\alpha_j I)q_j - \\beta_{j\\text{-}1} q_j\n\nIn other words, if (q_{j\\text{-}1}, \\beta_j, q_j) are known, then (\\alpha_j, \\beta_{j+1}, q_{j+1}) are completely determined. This fact is fantastic from a computational point of view: no explicit call to the QR algorithm necessary3!\nNote that a symmetric tridiagonal matrix is fully characterized by its diagonal and subdiagonal terms, which requires just O(n) space. If we assume that v \\mapsto Av \\sim O(n), then the above procedure clearly takes at most O(n^2) time, since there are most n such vectors \\{q_i\\}_{i=1}^n to generate!\nMoreover, if we only need the eigen-values \\Lambda(A) ( and not their vectors U), then we may execute the recurrence keeping at most three vectors \\{q_{j-1}, q_{j}, q_{j+1}\\} in memory at any given time. Since each of these is O(n) is size, the claim of O(n) space is justified!\nThe description above is essentially the proof of the following Theorem:\n\nTheorem 1 (Parlett 1994, Simon 1984) Given a symmetric rank-r matrix A \\in \\mathbb{R}^{n \\times n} whose operator x \\mapsto A x requires O(\\eta) time and O(\\nu) space, the Lanczos iteration computes \\Lambda(A) = \\{ \\lambda_1, \\lambda_2, \\dots, \\lambda_r \\} in O(\\max\\{\\eta, n\\}\\cdot r) time and O(\\max\\{\\nu, n\\}) space, when computation is done in exact arithmetic",
+    "text": "Quadratic time and linear space? How?\nUnless you know the tricks, its not obvious at all decompositions above takes just O(n^2) time and O(n) space to obtain. So… how does the complexity argument play out?\n\nThe cubic bound\nFirst, its important to establish that computing the eigen-decomposition A = U \\Lambda U^T for general symmetric A \\in \\mathbb{R}^{n \\times n} is essentially bounded by \\Theta(n^\\omega) time and \\Theta(n^2) space, where \\omega \\approx 2.37\\dots is the matrix-multiplication constant. Conceptually, if we exclude the Strassen-model for matrix multiplication (since it is not practical anyways), this translates to an effective \\Omega(n^3) time bound.\n\n\nSolving the right problem\nIn some applications, the eigenvectors are not needed all at once (or at all, even). One of the main draws to the Lanczos method is its efficiency: if one can perform v \\mapsto Av quickly—say, in \\approx O(n) time—then the Lanczos method can construct \\Lambda(A) in just O(n^2) time and O(n) space! Moreover, entire method is matrix free as the only input to the algorithm is a (fast) matrix-vector product v \\mapsto Av: one need not store A explicitly to do this for many special types of linear operators.\nIf you squint hard enough, you can deduce that since A Q_j = Q_j T_j + \\beta_{j+1} q_{j+1} e_{j}^T, every symmetric A \\in \\mathbb{R}^{n \\times n} expanded this way admits a three-term recurrence: \n\\begin{align*}\nA q_j &= \\beta_{j\\text{-}1} q_{j\\text{-}1} + \\alpha_j q_j + \\beta_j q_{j+1} \\\\\n\\Leftrightarrow \\beta_{j} \\, q_{j+1} &= A q_j - \\alpha_j \\, q_j - \\beta_{j\\text{-}1} \\, q_{j\\text{-}1}  \n\\end{align*}\n\nThe equation above is a variable-coefficient second-order linear difference equation, and it is known such equations have unique solutions: \n\\alpha_j = q_j^T A q_j, \\;\\; \\beta_j = \\lVert r_j \\rVert_2, \\;\\; q_{j+1} = r_j / \\beta_j\n\n\n\\text{where  } r_j = (A - \\alpha_j I)q_j - \\beta_{j\\text{-}1} q_j\n\nIn other words, if (q_{j\\text{-}1}, \\beta_j, q_j) are known, then (\\alpha_j, \\beta_{j+1}, q_{j+1}) are completely determined. This fact is fantastic from a computational point of view: no explicit call to the QR algorithm necessary3!\nNote that a symmetric tridiagonal matrix is fully characterized by its diagonal and subdiagonal terms, which requires just O(n) space. If we assume that v \\mapsto Av \\sim O(n), then the above procedure clearly takes at most O(n^2) time, since there are most n such vectors \\{q_i\\}_{i=1}^n to generate!\nMoreover, if we only need the eigen-values \\Lambda(A) ( and not their vectors U), then we may execute the recurrence keeping at most three vectors \\{q_{j-1}, q_{j}, q_{j+1}\\} in memory at any given time. Since each of these is O(n) is size, the claim of O(n) space is justified!\nThe description above is essentially the proof of the following Theorem:\n\nTheorem 1 (Parlett 1994, Simon 1984) Given a symmetric rank-r matrix A \\in \\mathbb{R}^{n \\times n} whose operator x \\mapsto A x requires O(\\eta) time and O(\\nu) space, the Lanczos iteration computes \\Lambda(A) = \\{ \\lambda_1, \\lambda_2, \\dots, \\lambda_r \\} in O(\\max\\{\\eta, n\\}\\cdot r) time and O(\\max\\{\\nu, n\\}) space, when computation is done in exact arithmetic",
     "crumbs": [
       "Theory",
       "Lanczos"

docs/src/basic/integration.qmd

@@ -11,5 +11,5 @@ More generally, the basic requirements for any operator `A` to be used with e.g.
 
 <!-- Note that in the matrix setting, `A.shape()` yields `(A.nrow(), A.ncol())`, so existing matrix implementations like [ndarray](https://numpy.org/doc/stable/reference/generated/numpy.ndarray.html)'s and [sparray](https://docs.scipy.org/doc/scipy/reference/sparse.html#sparse-array-classes)'s natively support the interface. -->
 
-These semantics extend to the C++ side as well via _C++20 concepts_---see the [C++ integration guide](advanced/cpp_integration.qmd). If you're using pybind11 and you want to extend `primate`'s Python API to work natively with linear operator implemented in C++, see the [pybind11 integration guide](advanced/pybind11_integration.qmd). 
+These semantics extend to the C++ side as well via _C++20 concepts_---see the [C++ integration guide](../advanced/cpp_integration.qmd). If you're using pybind11 and you want to extend `primate`'s Python API to work natively with linear operator implemented in C++, see the [pybind11 integration guide](../advanced/pybind11_integration.qmd). 
 

docs/src/theory/lanczos.qmd

@@ -2,9 +2,11 @@
 title: "The Lanczos method"
 editor: 
     rendor-on-save: true
+bibliography: /references.bib
 ---
 
-In 1950, Cornelius Lanczos proposed the *method of minimized iterations*---now known as the *Lanczos method*---enabling the world to **iteratively** expose the spectrum of a symmetric $A$ via *tridiagonalization*. Despite its age, the Lanczos iteration remains the standard[^1] algorithm both for computing eigensets and for solving linear systems in the large-scale regime. 
+In 1950, Cornelius Lanczos proposed the *method of minimized iterations*---now known as the *Lanczos method*---enabling the world to **iteratively** expose the spectra of symmetric operators $A \in \mathbb{S}_n$ via *tridiagonalization*. Despite its age, the Lanczos iteration remains the standard algorithm[^1] both for computing eigensets and for solving linear systems in the large-scale regime. 
+Having intrinsic connections to the conjugate gradient method, the theory of orthogonal polynomials, and Gaussian quadrature, it's fair to say the Lanczos method is one of the most important numerical methods of all time. Indeed, _Krylov methods_ listed among the top 10 most influential algorithms of the 20th century in an [IEEE guest editorial](https://www.computer.org/csdl/magazine/cs/2000/01/c1022/13rRUxBJhBm).
 
 <!-- In fact, a [an IEEE guest editorial](https://www.computer.org/csdl/magazine/cs/2000/01/c1022/13rRUxBJhBm) elected Krylov Subspace Iteration Methods among the top 10 most influential algorithms of the 20th century---on essentially equal footing with other well known algorithms, such as e.g. the _Fast Fourier Transform_, the _Simplex algorithm_, and _QuickSort_. 
 
@@ -15,6 +17,7 @@ $$\mathcal{K}_j(A, v) \triangleq \mathrm{span}\{ v, Av, A^2 v, \dots, A^{j-1}v \
 <!-- Unfortunately, a naive implementation of the Lanczos method can be practically useless for computing eigenvalues accurately:  --> 
 <!-- When Cornelius first published the method,  -->
 
+`primate` offers a matrix-free  
 In this post, I'll describe the basic theory behind the Lanczos method, with a focus on its time and space complexity. In general, the Lanczos method is intrinsically connected to other mathematical constructions, e.g. orthogonal polynomials, gaussian quadrature, conjugate gradient---none of these will be discussed for now. 
 
 ## Lanczos on a bumper sticker
@@ -51,16 +54,18 @@ That's pretty fortunate, because the eigen-sets of $T_j$ can be obtained in just
 
 To quote the lucid [Lanczos introduction from Parlett](https://apps.dtic.mil/sti/tr/pdf/ADA289614.pdf), _could anything be more simple?_
 
-## Surpassing the cubic bound
+## Quadratic time and linear space? How?
 
-Computing the eigen-decomposition $A = U \Lambda U^T$ for general symmetric $A \in \mathbb{R}^{n \times n}$ is essentially bounded by $\Theta(n^\omega)$ time and $\Theta(n^2)$ space, where $\omega \approx 2.37\dots$ is the matrix-multiplication constant. Conceptually, if we exclude the Strassen-model for matrix multiplication (since it is not practical anyways), this translates to an effective $\Omega(n^3)$ time bound. Not great!
+Unless you know the tricks, its not obvious at all decompositions above takes just $O(n^2)$ time and $O(n)$ space to obtain. So... how does the complexity argument play out?
 
-One of the main draws to the Lanczos method is its efficiency: if one can perform $v \mapsto Av$ quickly---say, in $\approx O(n)$ time---then the Lanczos method can construct $\Lambda(A)$ in _just_ [$O(n^2)$ time]{style="color: red;"} and [$O(n)$ space]{style="color: red;"}! 
-Moreover, entire method is *matrix free*     as the only input to the algorithm is a (fast) matrix-vector product $v \mapsto Av$: one need not store $A$ explicitly to do this for many special types of linear operators. 
+### The cubic bound 
 
-## Quadratic time and linear space? How?
+First, its important to establish that computing the eigen-decomposition $A = U \Lambda U^T$ for general symmetric $A \in \mathbb{R}^{n \times n}$ is essentially bounded by $\Theta(n^\omega)$ time and $\Theta(n^2)$ space, where $\omega \approx 2.37\dots$ is the matrix-multiplication constant. Conceptually, if we exclude the Strassen-model for matrix multiplication (since it is not practical anyways), this translates to an effective $\Omega(n^3)$ time bound. 
 
-Unless you know the tricks, its not obvious at all decompositions above takes just $O(n^2)$ time and $O(n)$ space to obtain. So... how does the complexity argument play out?
+### Solving the right problem
+
+In some applications, the eigenvectors are not needed all at once (or at all, even). One of the main draws to the Lanczos method is its efficiency: if one can perform $v \mapsto Av$ quickly---say, in $\approx O(n)$ time---then the Lanczos method can construct $\Lambda(A)$ in _just_ [$O(n^2)$ time]{style="color: red;"} and [$O(n)$ space]{style="color: red;"}! 
+Moreover, entire method is *matrix free* as the only input to the algorithm is a (fast) matrix-vector product $v \mapsto Av$: one need not store $A$ explicitly to do this for many special types of linear operators. 
 
 If you squint hard enough, you can deduce that since $A Q_j = Q_j T_j + \beta_{j+1} q_{j+1} e_{j}^T$, every symmetric $A \in \mathbb{R}^{n \times n}$ expanded this way admits a *three-term recurrence*: 
 $$