Add documentation for truncated decomposition keyword arguments

.github/workflows/PR_comment.yml

@@ -0,0 +1,17 @@
+name: Docs preview comment
+on:
+  pull_request:
+    types: [labeled]
+
+permissions:
+  pull-requests: write
+jobs:
+  pr_comment:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Create PR comment
+        if: github.event_name == 'pull_request' && github.repository == github.event.pull_request.head.repo.full_name && github.event.label.name == 'documentation'
+        uses: thollander/actions-comment-pull-request@v3
+        with:
+          message: 'After the build completes, the updated documentation will be available [here](https://quantumkithub.github.io/MatrixAlgebraKit.jl/previews/PR${{ github.event.number }}/)'
+          comment-tag: 'preview-doc'

docs/make.jl

@@ -52,6 +52,7 @@ makedocs(;
             "user_interface/compositions.md",
             "user_interface/decompositions.md",
             "user_interface/truncations.md",
+            "user_interface/properties.md",
             "user_interface/matrix_functions.md",
         ],
         "Developer Interface" => "dev_interface.md",

docs/src/index.md

@@ -26,6 +26,7 @@ These operations typically follow some common skeleton, and here we go into a li
 
 ```@contents
 Pages = ["user_interface/compositions.md", "user_interface/decompositions.md",
-         "user_interface/truncations.md", "user_interface/matrix_functions.md"]
+         "user_interface/truncations.md", "user_interface/properties.md",
+         "user_interface/matrix_functions.md"]
 Depth = 2
 ```

docs/src/user_interface/decompositions.md

@@ -55,6 +55,7 @@ Not all matrices can be diagonalized, and some real matrices can only be diagona
 In particular, the resulting decomposition can only guaranteed to be real for real symmetric inputs `A`.
 Therefore, we provide `eig_` and `eigh_` variants, where `eig` always results in complex-valued `V` and `D`, while `eigh` requires symmetric inputs but retains the scalartype of the input.
 
+The full set of eigenvalues and eigenvectors can be computed using the [`eig_full`](@ref) and [`eigh_full`](@ref) functions.
 If only the eigenvalues are required, the [`eig_vals`](@ref) and [`eigh_vals`](@ref) functions can be used.
 These functions return the diagonal elements of `D` in a vector.
 
@@ -99,7 +100,7 @@ Filter = t -> t isa Type && t <: MatrixAlgebraKit.LAPACK_EigAlgorithm
 
 The [Schur decomposition](https://en.wikipedia.org/wiki/Schur_decomposition) transforms a complex square matrix `A` into a product `Q * T * Qᴴ`, where `Q` is unitary and `T` is upper triangular.
 It rewrites an arbitrary complex square matrix as unitarily similar to an upper triangular matrix whose diagonal elements are the eigenvalues of `A`.
-For real matrices, the same decomposition can be achieved with `T` being quasi-upper triangular, ie triangular with blocks of size `(1, 1)` and `(2, 2)` on the diagonal.
+For real matrices, the same decomposition can be achieved in real arithmetic by allowing `T` to be quasi-upper triangular, i.e. triangular with blocks of size `(1, 1)` and `(2, 2)` on the diagonal.
 
 This decomposition is also useful for computing the eigenvalues of a matrix, which is exposed through the [`schur_vals`](@ref) function.
 
@@ -117,12 +118,12 @@ Filter = t -> t isa Type && t <: MatrixAlgebraKit.LAPACK_EigAlgorithm
 
 ## Singular Value Decomposition
 
-The [Singular Value Decomposition](https://en.wikipedia.org/wiki/Singular_value_decomposition) transforms a matrix `A` into a product `U * Σ * Vᴴ`, where `U` and `V` are orthogonal, and `Σ` is diagonal, real and non-negative.
-For a square matrix `A`, both `U` and `V` are unitary, and if the singular values are distinct, the decomposition is unique.
+The [Singular Value Decomposition](https://en.wikipedia.org/wiki/Singular_value_decomposition) transforms a matrix `A` into a product `U * Σ * Vᴴ`, where `U` and `Vᴴ` are unitary, and `Σ` is diagonal, real and non-negative.
+For a square matrix `A`, both `U` and `Vᴴ` are unitary, and if the singular values are distinct, the decomposition is unique.
 
 For rectangular matrices `A` of size `(m, n)`, there are two modes of operation, [`svd_full`](@ref) and [`svd_compact`](@ref).
-The former ensures that the resulting `U`, and `V` remain square unitary matrices, of size `(m, m)` and `(n, n)`, with rectangular `Σ` of size `(m, n)`.
-The latter creates an isometric `U` of size `(m, min(m, n))`, and `V` of size `(n, min(m, n))`, with a square `Σ` of size `(min(m, n), min(m, n))`.
+The former ensures that the resulting `U`, and `Vᴴ` remain square unitary matrices, of size `(m, m)` and `(n, n)`, with rectangular `Σ` of size `(m, n)`.
+The latter creates an isometric `U` of size `(m, min(m, n))`, and `V = (Vᴴ)'` of size `(n, min(m, n))`, with a square `Σ` of size `(min(m, n), min(m, n))`.
 
 It is also possible to compute the singular values only, using the [`svd_vals`](@ref) function.
 This then returns a vector of the values on the diagonal of `Σ`.
@@ -147,21 +148,23 @@ Filter = t -> t isa Type && t <: MatrixAlgebraKit.LAPACK_SVDAlgorithm
 
 The [Polar Decomposition](https://en.wikipedia.org/wiki/Polar_decomposition) of a matrix `A` is a factorization `A = W * P`, where `W` is unitary and `P` is positive semi-definite.
 If `A` is invertible (and therefore square), the polar decomposition always exists and is unique.
-For non-square matrices, the polar decomposition is not unique, but `P` is.
-In particular, the polar decomposition is unique if `A` is full rank.
+For non-square matrices `A` of size `(m, n)`, the decomposition `A = W * P` with `P` positive semi-definite of size `(n, n)` and `W` isometric of size `(m, n)` exists only if `m >= n`, and is unique if `A` and thus `P` is full rank.
+For `m <= n`, we can analoguously decompose `A` as `A = P * Wᴴ` with `P` positive semi-definite of size `(m, m)` and `Wᴴ` of size `(m, n)` such that `W = (Wᴴ)'` is isometric. Only in the case `m = n` do both decompositions exist.
 
-This decomposition can be computed for both sides, resulting in the [`left_polar`](@ref) and [`right_polar`](@ref) functions.
+The decompositions `A = W * P` or `A = P * Wᴴ` can be computed with the [`left_polar`](@ref) and [`right_polar`](@ref) functions, respectively.
 
 ```@docs; canonical=false
 left_polar
 right_polar
 ```
 
-These functions are implemented by first computing a singular value decomposition, and then constructing the polar decomposition from the singular values and vectors.
-Therefore, the relevant LAPACK-based implementation is the one for the SVD:
+These functions can be implemented by first computing a singular value decomposition, and then constructing the polar decomposition from the singular values and vectors. Alternatively, the polar decomposition can be computed using an 
+iterative method based on Newton's method, that can be more efficient for large matrices, especially if they are
+close to being isometric already.
 
 ```@docs; canonical=false
 PolarViaSVD
+PolarNewton
 ```
 
 ## Orthogonal Subspaces
@@ -179,7 +182,7 @@ right_orth
 ## Null Spaces
 
 Similarly, it can be convenient to obtain an orthogonal basis for the kernel or cokernel of a matrix.
-These are the compliments of the image and coimage, and can be computed using the [`left_null`](@ref) and [`right_null`](@ref) functions.
+These are the compliments of the coimage and image, respectively, and can be computed using the [`left_null`](@ref) and [`right_null`](@ref) functions.
 Again, this is typically implemented through a combination of the decompositions mentioned above, and serves as a convenient interface to these operations.
 
 ```@docs; canonical=false

docs/src/user_interface/properties.md

@@ -0,0 +1,23 @@
+```@meta
+CurrentModule = MatrixAlgebraKit
+CollapsedDocStrings = true
+```
+
+# Matrix Properties
+
+MatrixAlgebraKit.jl provides a number of methods to check various properties of matrices.
+
+```@docs; canonical=false
+isisometric
+isunitary
+ishermitian
+isantihermitian
+```
+
+Furthermore, there are also methods to project a matrix onto the nearest matrix (in 2-norm distance) with a given property.
+
+```@docs; canonical=false
+project_isometric
+project_hermitian
+project_antihermitian
+```

docs/src/user_interface/truncations.md

@@ -5,7 +5,112 @@ CollapsedDocStrings = true
 
 # Truncations
 
-Currently, truncations are supported through the following different methods:
+Truncation strategies allow you to control which eigenvalues or singular values to keep when computing partial or truncated decompositions. These strategies are used in the functions [`eigh_trunc`](@ref), [`eig_trunc`](@ref), and [`svd_trunc`](@ref) to reduce the size of the decomposition while retaining the most important information.
+
+## Using Truncations in Decompositions
+
+Truncation strategies can be used with truncated decomposition functions in two ways, as illustrated below.
+For concreteness, we use the following matrix as an example:
+
+```jldoctest truncations
+using MatrixAlgebraKit
+using MatrixAlgebraKit: diagview
+
+A = [2 1 0; 1 3 1; 0 1 4];
+D, V = eigh_full(A);
+
+diagview(D) ≈ [3 - √3, 3, 3 + √3]
+
+# output
+
+true
+```
+
+### 1. Using the `trunc` keyword with a `NamedTuple`
+
+The simplest approach is to pass a `NamedTuple` with the truncation parameters.
+For example, keeping only the largest 2 eigenvalues:
+
+```jldoctest truncations
+Dtrunc, Vtrunc = eigh_trunc(A; trunc = (maxrank = 2,));
+size(Dtrunc, 1) <= 2
+
+# output
+
+true
+```
+
+Note however that there are no guarantees on the order of the output values:
+
+```jldoctest truncations
+diagview(Dtrunc) ≈ diagview(D)[[3, 2]]
+
+# output
+
+true
+```
+
+You can also use tolerance-based truncation or combine multiple criteria:
+
+```jldoctest truncations
+Dtrunc, Vtrunc = eigh_trunc(A; trunc = (atol = 2.9,));
+all(>(2.9), diagview(Dtrunc))
+
+# output
+
+true
+```
+
+```jldoctest truncations
+Dtrunc, Vtrunc = eigh_trunc(A; trunc = (maxrank = 2, atol = 2.9));
+size(Dtrunc, 1) <= 2 && all(>(2.9), diagview(Dtrunc))
+
+# output
+true
+```
+
+In general, the keyword arguments that are supported can be found in the `TruncationStrategy` docstring:
+
+```@docs; canonical = false
+TruncationStrategy
+```
+
+
+### 2. Using explicit `TruncationStrategy` objects
+
+For more control, you can construct [`TruncationStrategy`](@ref) objects directly.
+This is also what the previous syntax will end up calling.
+
+```jldoctest truncations
+Dtrunc, Vtrunc = eigh_trunc(A; trunc = truncrank(2))
+size(Dtrunc, 1) <= 2
+
+# output
+
+true
+```
+
+```jldoctest truncations
+Dtrunc, Vtrunc = eigh_trunc(A; trunc = truncrank(2) & trunctol(; atol = 2.9))
+size(Dtrunc, 1) <= 2 && all(>(2.9), diagview(Dtrunc))
+
+# output
+true
+```
+
+## Truncation with SVD vs Eigenvalue Decompositions
+
+When using truncations with different decomposition types, keep in mind:
+
+- **`svd_trunc`**: Singular values are always real and non-negative, sorted in descending order. Truncation by value typically keeps the largest singular values.
+
+- **`eigh_trunc`**: Eigenvalues are real but can be negative for symmetric matrices. By default, `truncrank` sorts by absolute value, so `truncrank(k)` keeps the `k` eigenvalues with largest magnitude (positive or negative).
+
+- **`eig_trunc`**: For general (non-symmetric) matrices, eigenvalues can be complex. Truncation by absolute value considers the complex magnitude.
+
+## Truncation Strategies
+
+MatrixAlgebraKit provides several built-in truncation strategies:
 
 ```@docs; canonical=false
 notrunc
@@ -15,11 +120,10 @@ truncfilter
 truncerror
 ```
 
-It is additionally possible to combine truncation strategies by making use of the `&` operator.
-For example, truncating to a maximal dimension `10`, and discarding all values below `1e-6` would be achieved by:
+Truncation strategies can be combined using the `&` operator to create intersection-based truncation.
+When strategies are combined, only the values that satisfy all conditions are kept.
 
 ```julia
-maxdim = 10
-atol = 1e-6
-combined_trunc = truncrank(maxdim) & trunctol(; atol)
+combined_trunc = truncrank(10) & trunctol(; atol = 1e-6);
 ```
+

src/interface/eig.jl

@@ -3,9 +3,10 @@
 # TODO: kwargs for sorting eigenvalues?
 
 docs_eig_note = """
-Note that [`eig_full`](@ref) and its variants do not assume additional structure on the input,
-and therefore will always return complex eigenvalues and eigenvectors. For the real
-eigenvalue decomposition of symmetric or hermitian operators, see [`eigh_full`](@ref).
+    Note that [`eig_full`](@ref) and its variants do not assume any symmetry structure on
+    the input matrices, and therefore will always return complex eigenvalues and eigenvectors
+    for reasons of type stability.  For the  eigenvalue decomposition of symmetric or hermitian
+    operators, see [`eigh_full`](@ref).
 """
 
 """
@@ -24,32 +25,53 @@ and the diagonal matrix `D` contains the associated eigenvalues.
     as it may not always be possible to use the provided `DV` as output.
 
 !!! note
-    $(docs_eig_note)
+$(docs_eig_note)
 
 See also [`eig_vals(!)`](@ref eig_vals) and [`eig_trunc(!)`](@ref eig_trunc).
 """
 @functiondef eig_full
 
 """
-    eig_trunc(A; kwargs...) -> D, V
+    eig_trunc(A; [trunc], kwargs...) -> D, V
     eig_trunc(A, alg::AbstractAlgorithm) -> D, V
-    eig_trunc!(A, [DV]; kwargs...) -> D, V
+    eig_trunc!(A, [DV]; [trunc], kwargs...) -> D, V
     eig_trunc!(A, [DV], alg::AbstractAlgorithm) -> D, V
 
 Compute a partial or truncated eigenvalue decomposition of the matrix `A`,
 such that `A * V ≈ V * D`, where the (possibly rectangular) matrix `V` contains 
 a subset of eigenvectors and the diagonal matrix `D` contains the associated eigenvalues,
 selected according to a truncation strategy.
 
+## Keyword arguments
+The behavior of this function is controlled by the following keyword arguments:
+
+- `trunc`: Specifies the truncation strategy. This can be:
+  - A `NamedTuple` with fields `atol`, `rtol`, and/or `maxrank`, which will be converted to
+    a [`TruncationStrategy`](@ref). For details on available truncation strategies, see
+    [Truncations](@ref).
+  - A `TruncationStrategy` object directly (e.g., `truncrank(10)`, `trunctol(atol=1e-6)`, or
+    combinations using `&`).
+  - `nothing` (default), which keeps all eigenvalues.
+
+- Other keyword arguments are passed to the algorithm selection procedure. If no explicit
+  `alg` is provided, these keywords are used to select and configure the algorithm through
+  [`MatrixAlgebraKit.select_algorithm`](@ref). The remaining keywords after algorithm
+  selection are passed to the algorithm constructor. See [`MatrixAlgebraKit.default_algorithm`](@ref)
+  for the default algorithm selection behavior.
+
+When `alg` is a [`TruncatedAlgorithm`](@ref), the `trunc` keyword cannot be specified as the
+truncation strategy is already embedded in the algorithm.
+
 !!! note
     The bang method `eig_trunc!` optionally accepts the output structure and
     possibly destroys the input matrix `A`. Always use the return value of the function
     as it may not always be possible to use the provided `DV` as output.
 
 !!! note
-    $(docs_eig_note)
+$(docs_eig_note)
 
-See also [`eig_full(!)`](@ref eig_full) and [`eig_vals(!)`](@ref eig_vals).
+See also [`eig_full(!)`](@ref eig_full), [`eig_vals(!)`](@ref eig_vals), and
+[Truncations](@ref) for more information on truncation strategies.
 """
 @functiondef eig_trunc
 
@@ -67,7 +89,7 @@ Compute the list of eigenvalues of `A`.
     as it may not always be possible to use the provided `D` as output.
 
 !!! note
-    $(docs_eig_note)
+$(docs_eig_note)
 
 See also [`eig_full(!)`](@ref eig_full) and [`eig_trunc(!)`](@ref eig_trunc).
 """