Merge branch 'main' into jh/truncerr

Author: lkdvos

.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);
 ```
+

ext/MatrixAlgebraKitChainRulesCoreExt.jl

@@ -14,7 +14,7 @@ MatrixAlgebraKit.iszerotangent(::AbstractZero) = true
 @non_differentiable MatrixAlgebraKit.select_algorithm(args...)
 @non_differentiable MatrixAlgebraKit.initialize_output(args...)
 @non_differentiable MatrixAlgebraKit.check_input(args...)
-@non_differentiable MatrixAlgebraKit.isisometry(args...)
+@non_differentiable MatrixAlgebraKit.isisometric(args...)
 @non_differentiable MatrixAlgebraKit.isunitary(args...)
 
 function ChainRulesCore.rrule(::typeof(copy_input), f, A)

src/MatrixAlgebraKit.jl

@@ -9,7 +9,7 @@ using LinearAlgebra: Diagonal, diag, diagind, isdiag
 using LinearAlgebra: UpperTriangular, LowerTriangular
 using LinearAlgebra: BlasFloat, BlasReal, BlasComplex, BlasInt
 
-export isisometry, isunitary, ishermitian, isantihermitian
+export isisometric, isunitary, ishermitian, isantihermitian
 
 export project_hermitian, project_antihermitian, project_isometric
 export project_hermitian!, project_antihermitian!, project_isometric!
@@ -65,6 +65,7 @@ export notrunc, truncrank, trunctol, truncerror, truncfilter
             :svd_pullback!, :svd_trunc_pullback!
         )
     )
+    eval(Expr(:public, :is_left_isometric, :is_right_isometric))
 end
 
 include("common/defaults.jl")

src/common/matrixproperties.jl

@@ -1,5 +1,5 @@
 """
-    isisometry(A; side=:left, isapprox_kwargs...) -> Bool
+    isisometric(A; side=:left, isapprox_kwargs...) -> Bool
 
 Test whether a linear map is an isometry, where the type of isometry is controlled by `kind`:
 
@@ -8,13 +8,14 @@ Test whether a linear map is an isometry, where the type of isometry is controll
 
 The `isapprox_kwargs` are passed on to `isapprox` to control the tolerances.
 
-New specializations should overload [`is_left_isometry`](@ref) and [`is_right_isometry`](@ref).
+New specializations should overload [`MatrixAlgebraKit.is_left_isometric`](@ref) and
+[`MatrixAlgebraKit.is_right_isometric`](@ref).
 
 See also [`isunitary`](@ref).
 """
-function isisometry(A; side::Symbol = :left, isapprox_kwargs...)
-    side === :left && return is_left_isometry(A; isapprox_kwargs...)
-    side === :right && return is_right_isometry(A; isapprox_kwargs...)
+function isisometric(A; side::Symbol = :left, isapprox_kwargs...)
+    side === :left && return is_left_isometric(A; isapprox_kwargs...)
+    side === :right && return is_right_isometric(A; isapprox_kwargs...)
 
     throw(ArgumentError(lazy"Invalid isometry side: $side"))
 end
@@ -25,48 +26,42 @@ end
 Test whether a linear map is unitary, i.e. `A * A' ≈ I ≈ A' * A`.
 The `isapprox_kwargs` are passed on to `isapprox` to control the tolerances.
 
-See also [`isisometry`](@ref).
+See also [`isisometric`](@ref).
 """
 function isunitary(A; isapprox_kwargs...)
-    return is_left_isometry(A; isapprox_kwargs...) &&
-        is_right_isometry(A; isapprox_kwargs...)
+    return is_left_isometric(A; isapprox_kwargs...) &&
+        is_right_isometric(A; isapprox_kwargs...)
 end
 function isunitary(A::AbstractMatrix; isapprox_kwargs...)
     size(A, 1) == size(A, 2) || return false
-    return is_left_isometry(A; isapprox_kwargs...)
+    return is_left_isometric(A; isapprox_kwargs...)
 end
 
 @doc """
-    is_left_isometry(A; isapprox_kwargs...) -> Bool
+    is_left_isometric(A; isapprox_kwargs...) -> Bool
 
-Test whether a linear map is a left isometry, i.e. `A' * A ≈ I`.
+Test whether a linear map is a (left) isometry, i.e. `A' * A ≈ I`.
 The `isapprox_kwargs` can be used to control the tolerances of the equality.
 
-See also [`isisometry`](@ref) and [`is_right_isometry`](@ref).
-""" is_left_isometry
+See also [`isisometric`](@ref) and [`MatrixAlgebraKit.is_right_isometric`](@ref).
+""" is_left_isometric
 
-function is_left_isometry(A::AbstractMatrix; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
+function is_left_isometric(A::AbstractMatrix; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
     P = A' * A
     nP = norm(P) # isapprox would use `rtol * max(norm(P), norm(I))`
     diagview(P) .-= 1
     return norm(P) <= max(atol, rtol * nP) # assume that the norm of I is `sqrt(n)`
 end
 
 @doc """
-    is_right_isometry(A; isapprox_kwargs...) -> Bool
+    is_right_isometric(A; isapprox_kwargs...) -> Bool
 
-Test whether a linear map is a right isometry, i.e. `A * A' ≈ I`.
+Test whether a linear map is a (right) isometry, i.e. `A * A' ≈ I`.
 The `isapprox_kwargs` can be used to control the tolerances of the equality.
 
-See also [`isisometry`](@ref) and [`is_left_isometry`](@ref).
-""" is_right_isometry
-
-function is_right_isometry(A::AbstractMatrix; atol::Real = 0, rtol::Real = defaulttol(A), norm = LinearAlgebra.norm)
-    P = A * A'
-    nP = norm(P) # isapprox would use `rtol * max(norm(P), norm(I))`
-    diagview(P) .-= 1
-    return norm(P) <= max(atol, rtol * nP) # assume that the norm of I is `sqrt(n)`
-end
+See also [`isisometric`](@ref) and [`MatrixAlgebraKit.is_left_isometric`](@ref).
+""" is_right_isometric
+is_right_isometric(A; kwargs...) = is_left_isometric(A'; kwargs...)
 
 """
     ishermitian(A; isapprox_kwargs...)