diff --git a/.github/workflows/PR_comment.yml b/.github/workflows/PR_comment.yml new file mode 100644 index 00000000..fd1f1d5f --- /dev/null +++ b/.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' diff --git a/docs/make.jl b/docs/make.jl index 3e37770b..c2ae1fcc 100644 --- a/docs/make.jl +++ b/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", diff --git a/docs/src/index.md b/docs/src/index.md index 20e3185a..7055ba3e 100644 --- a/docs/src/index.md +++ b/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 ``` diff --git a/docs/src/user_interface/decompositions.md b/docs/src/user_interface/decompositions.md index 5a25d3ce..b6d0e94e 100644 --- a/docs/src/user_interface/decompositions.md +++ b/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 diff --git a/docs/src/user_interface/properties.md b/docs/src/user_interface/properties.md new file mode 100644 index 00000000..badf4281 --- /dev/null +++ b/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 +``` \ No newline at end of file diff --git a/docs/src/user_interface/truncations.md b/docs/src/user_interface/truncations.md index 40799890..e6038d32 100644 --- a/docs/src/user_interface/truncations.md +++ b/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); ``` + diff --git a/src/interface/eig.jl b/src/interface/eig.jl index 5a74f79f..62dc713e 100644 --- a/src/interface/eig.jl +++ b/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,16 +25,16 @@ 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`, @@ -41,15 +42,36 @@ such that `A * V ≈ V * D`, where the (possibly rectangular) matrix `V` contain 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). """ diff --git a/src/interface/eigh.jl b/src/interface/eigh.jl index c51b2c01..f1fe1936 100644 --- a/src/interface/eigh.jl +++ b/src/interface/eigh.jl @@ -3,9 +3,12 @@ # TODO: kwargs for sorting eigenvalues? docs_eigh_note = """ -Note that [`eigh_full`](@ref) and its variants assume additional structure on the input, -and therefore will retain the `eltype` of the input for the eigenvalues and eigenvectors. -For generic eigenvalue decompositions, see [`eig_full`](@ref). + Note that [`eigh_full`](@ref) and its variants assume that the input matrix is hermitian, + or thus symmetric if the input is real. The resulting algorithms exploit this structure, + and return eigenvalues that are always real, and eigenvectors that are orthogonal and have + the same `eltype` as the input matrix. If the input matrix does not have this structure, + the generic eigenvalue decomposition provided by [`eig_full`](@ref) and its variants + should be used instead. """ """ @@ -24,22 +27,42 @@ and the real diagonal matrix `D` contains the associated eigenvalues. as it may not always be possible to use the provided `DV` as output. !!! note - $(docs_eigh_note) +$(docs_eigh_note) See also [`eigh_vals(!)`](@ref eigh_vals) and [`eigh_trunc(!)`](@ref eigh_trunc). """ @functiondef eigh_full """ - eigh_trunc(A; kwargs...) -> D, V + eigh_trunc(A; [trunc], kwargs...) -> D, V eigh_trunc(A, alg::AbstractAlgorithm) -> D, V - eigh_trunc!(A, [DV]; kwargs...) -> D, V + eigh_trunc!(A, [DV]; [trunc], kwargs...) -> D, V eigh_trunc!(A, [DV], alg::AbstractAlgorithm) -> D, V Compute a partial or truncated eigenvalue decomposition of the symmetric or hermitian matrix `A`, such that `A * V ≈ V * D`, where the isometric matrix `V` contains a subset of the orthogonal eigenvectors and the real diagonal matrix `D` contains the associated eigenvalues, -selected according to a truncation strategy. +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 `eigh_trunc!` optionally accepts the output structure and @@ -47,9 +70,10 @@ selected according to a truncation strategy. as it may not always be possible to use the provided `DV` as output. !!! note - $(docs_eigh_note) +$(docs_eigh_note) -See also [`eigh_full(!)`](@ref eigh_full) and [`eigh_vals(!)`](@ref eigh_vals). +See also [`eigh_full(!)`](@ref eigh_full), [`eigh_vals(!)`](@ref eigh_vals), and +[Truncations](@ref) for more information on truncation strategies. """ @functiondef eigh_trunc @@ -67,7 +91,7 @@ Compute the list of (real) eigenvalues of the symmetric or hermitian matrix `A`. as it may not always be possible to use the provided `DV` as output. !!! note - $(docs_eigh_note) +$(docs_eigh_note) See also [`eigh_full(!)`](@ref eigh_full) and [`eigh_trunc(!)`](@ref eigh_trunc). """ diff --git a/src/interface/gen_eig.jl b/src/interface/gen_eig.jl index 5b7cfc88..92cd3477 100644 --- a/src/interface/gen_eig.jl +++ b/src/interface/gen_eig.jl @@ -3,9 +3,11 @@ # TODO: kwargs for sorting eigenvalues? docs_gen_eig_note = """ -Note that [`gen_eig_full`](@ref) and its variants do not assume additional structure on the inputs, -and therefore will always return complex eigenvalues and eigenvectors. For the real -generalized eigenvalue decomposition is not yet supported. + Note that [`gen_eig_full`](@ref) and its variants do not assume additional structure + on the inputs, and therefore will always return complex eigenvalues and eigenvectors. + Real eigenvalues can be expected when both input matrices are hermitian and one of them + is positive definite, but specialized methods that exploit this structure are not yet + implemented or supported. """ # TODO: do we need "full"? @@ -26,7 +28,7 @@ and the diagonal matrix `W` contains the associated generalized eigenvalues. possible to use the provided `WV` as output. !!! note - $(docs_gen_eig_note) +$(docs_gen_eig_note) See also [`gen_eig_vals(!)`](@ref eig_vals). """ @@ -47,7 +49,7 @@ Compute the list of generalized eigenvalues of `A` and `B`. provided `W` as output. !!! note - $(docs_gen_eig_note) +$(docs_gen_eig_note) See also [`gen_eig_full(!)`](@ref gen_eig_full). """ diff --git a/src/interface/svd.jl b/src/interface/svd.jl index ee6cd9af..e36b0462 100644 --- a/src/interface/svd.jl +++ b/src/interface/svd.jl @@ -41,11 +41,10 @@ See also [`svd_full(!)`](@ref svd_full), [`svd_vals(!)`](@ref svd_vals) and """ @functiondef svd_compact -# TODO: decide if we should have `svd_trunc!!` instead """ - svd_trunc(A; kwargs...) -> U, S, Vᴴ + svd_trunc(A; [trunc], kwargs...) -> U, S, Vᴴ svd_trunc(A, alg::AbstractAlgorithm) -> U, S, Vᴴ - svd_trunc!(A, [USVᴴ]; kwargs...) -> U, S, Vᴴ + svd_trunc!(A, [USVᴴ]; [trunc], kwargs...) -> U, S, Vᴴ svd_trunc!(A, [USVᴴ], alg::AbstractAlgorithm) -> U, S, Vᴴ Compute a partial or truncated singular value decomposition (SVD) of `A`, such that @@ -53,14 +52,34 @@ Compute a partial or truncated singular value decomposition (SVD) of `A`, such t `(m, k)`, whereas `Vᴴ` is a matrix of size `(k, n)` with orthonormal rows and `S` is a square diagonal matrix of size `(k, k)`, with `k` is set by the 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 singular values. + +- 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 `svd_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 `USVᴴ` as output. - -See also [`svd_full(!)`](@ref svd_full), [`svd_compact(!)`](@ref svd_compact) and -[`svd_vals(!)`](@ref svd_vals). +See also [`svd_full(!)`](@ref svd_full), [`svd_compact(!)`](@ref svd_compact), +[`svd_vals(!)`](@ref svd_vals), and [Truncations](@ref) for more information on +truncation strategies. """ @functiondef svd_trunc diff --git a/src/interface/truncation.jl b/src/interface/truncation.jl index 9504b48f..10ace743 100644 --- a/src/interface/truncation.jl +++ b/src/interface/truncation.jl @@ -4,26 +4,36 @@ Select a truncation strategy based on the provided keyword arguments. ## Keyword arguments -- `atol=nothing` : Absolute tolerance for the truncation -- `rtol=nothing` : Relative tolerance for the truncation -- `maxrank=nothing` : Maximal rank for the truncation -""" -function TruncationStrategy(; atol = nothing, rtol = nothing, maxrank = nothing) - if isnothing(maxrank) && isnothing(atol) && isnothing(rtol) - return NoTruncation() - elseif isnothing(maxrank) +The following keyword arguments are all optional, and their default value (`nothing`) +will be ignored. It is also allowed to combine multiple of these, in which case the kept +values will consist of the intersection of the different truncated strategies. + +- `atol::Real` : Absolute tolerance for the truncation +- `rtol::Real` : Relative tolerance for the truncation +- `maxrank::Real` : Maximal rank for the truncation +- `maxerror::Real` : Maximal truncation error. +- `filter` : Custom filter to select truncated values. +""" +function TruncationStrategy(; + atol::Union{Real, Nothing} = nothing, + rtol::Union{Real, Nothing} = nothing, + maxrank::Union{Real, Nothing} = nothing, + maxerror::Union{Real, Nothing} = nothing, + filter = nothing + ) + strategy = notrunc() + + if !isnothing(atol) || !isnothing(rtol) atol = @something atol 0 rtol = @something rtol 0 - return trunctol(; atol, rtol) - else - if isnothing(atol) && isnothing(rtol) - return truncrank(maxrank) - else - atol = @something atol 0 - rtol = @something rtol 0 - return truncrank(maxrank) & trunctol(; atol, rtol) - end + strategy &= trunctol(; atol, rtol) end + + isnothing(maxrank) || (strategy &= truncrank(maxrank)) + isnothing(maxerror) || (strategy &= truncerror(; atol = maxerror)) + isnothing(filter) || (strategy &= truncfilter(filter)) + + return strategy end """ @@ -151,6 +161,8 @@ end function Base.:&(trunc1::TruncationStrategy, trunc2::TruncationStrategy) return TruncationIntersection((trunc1, trunc2)) end + +# flatten components function Base.:&(trunc1::TruncationIntersection, trunc2::TruncationIntersection) return TruncationIntersection((trunc1.components..., trunc2.components...)) end @@ -160,3 +172,12 @@ end function Base.:&(trunc1::TruncationStrategy, trunc2::TruncationIntersection) return TruncationIntersection((trunc1, trunc2.components...)) end + +# drop notrunc +Base.:&(::NoTruncation, trunc::TruncationStrategy) = trunc +Base.:&(trunc::TruncationStrategy, ::NoTruncation) = trunc +Base.:&(::NoTruncation, ::NoTruncation) = notrunc() + +# disambiguate +Base.:&(::NoTruncation, trunc::TruncationIntersection) = trunc +Base.:&(trunc::TruncationIntersection, ::NoTruncation) = trunc diff --git a/test/truncate.jl b/test/truncate.jl index 94a788e3..f2d24a7a 100644 --- a/test/truncate.jl +++ b/test/truncate.jl @@ -28,7 +28,7 @@ using MatrixAlgebraKit: NoTruncation, TruncationIntersection, TruncationByOrder, trunc = @constinferred TruncationStrategy(; atol, rtol, maxrank) @test trunc isa TruncationIntersection - @test trunc == truncrank(maxrank) & trunctol(; atol, rtol) + @test trunc == trunctol(; atol, rtol) & truncrank(maxrank) values = [1, 0.9, 0.5, -0.3, 0.01] @test values[@constinferred(findtruncated(values, truncrank(2)))] == values[1:2]