diff --git a/Project.toml b/Project.toml index a4e1c820..a7851a45 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MatrixAlgebraKit" uuid = "6c742aac-3347-4629-af66-fc926824e5e4" authors = ["Jutho and contributors"] -version = "0.5.0" +version = "0.6.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/docs/make.jl b/docs/make.jl index c2ae1fcc..4584ff1c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -57,6 +57,7 @@ makedocs(; ], "Developer Interface" => "dev_interface.md", "Library" => "library.md", + "Changelog" => "changelog.md", ], checkdocs = :exports, doctest = true, diff --git a/docs/src/changelog.md b/docs/src/changelog.md new file mode 100644 index 00000000..3eaab026 --- /dev/null +++ b/docs/src/changelog.md @@ -0,0 +1,64 @@ +# Changelog + +All notable changes to this project will be documented in this file. + +The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), +and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). + +## Guidelines for updating this changelog + +When making changes to this project, please update the "Unreleased" section with your changes under the appropriate category: + +- **Added** for new features. +- **Changed** for changes in existing functionality. +- **Deprecated** for soon-to-be removed features. +- **Removed** for now removed features. +- **Fixed** for any bug fixes. +- **Security** in case of vulnerabilities. + +When releasing a new version, move the "Unreleased" changes to a new version section with the release date. + +[Unreleased]: https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/compare/v0.6.0...HEAD +[0.6.0]: https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/releases/tag/v0.6.0 + +## [Unreleased] + +### Added + +### Changed + +### Deprecated + +### Removed + +### Fixed + +### Security + +## [0.6.0] - 2025-11-14 + +### Added +- New `project_isometric` function for projecting matrices onto isometric manifold ([#67](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/67)) +- New `PolarNewton` algorithm for polar decomposition ([#67](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/67)) +- New matrix property functions: `ishermitian`, `isantihermitian`, `hermitianpart!`, `hermitianpart`, `antihermitianpart!`, and `antihermitianpart` ([#64](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/64)) +- Support for `BigFloat` via new `GenericLinearAlgebra` extension ([#87](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/87)) +- Mooncake reverse-mode AD rules ([#85](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/85)) +- GPU support for image and null space computations ([#82](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/82)) +- GPU support for polar decomposition ([#83](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/83)) +- GPU support for new projection operations ([#81](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/81)) +- Output truncation error for truncated decompositions ([#75](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/75)) +- Documentation for truncated decomposition keyword arguments ([#71](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/71)) +- Default algorithm implementations for GPU wrapper array types ([#49](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/49)) + +### Changed + +- Made `gaugefix!` optional ([#95](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/95)) +- Renamed `isisometry` to `isisometric` for consistency with `project_isometric` ([#73](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/73)) +- Refactored `left_orth`, `right_orth`, `left_null` and `right_null` interface ([#79](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/79)) +- Improved GPU support for SVD operations ([#80](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/80)) +- Loosened strictness on hermitian checks ([#78](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/78)) +- Updated pullback tolerances ([#92](https://github.com/QuantumKitHub/MatrixAlgebraKit.jl/pull/92)) + +### Removed + +### Fixed diff --git a/docs/src/user_interface/decompositions.md b/docs/src/user_interface/decompositions.md index 60959cc4..fe28fcf9 100644 --- a/docs/src/user_interface/decompositions.md +++ b/docs/src/user_interface/decompositions.md @@ -413,21 +413,21 @@ This preserves the decomposition `A = U * Σ * Vᴴ` while fixing the gauge. ### Disabling Gauge Fixing Gauge fixing is enabled by default for all eigenvalue and singular value decompositions. -If you prefer to obtain the raw results from the underlying computational routines without gauge fixing, you can disable it using the `gaugefix` keyword argument: +If you prefer to obtain the raw results from the underlying computational routines without gauge fixing, you can disable it using the `fixgauge` keyword argument: ```julia # With gauge fixing (default) D, V = eigh_full(A) # Without gauge fixing -D, V = eigh_full(A; gaugefix = false) +D, V = eigh_full(A; fixgauge = false) ``` The same keyword is available for `eig_full`, `eig_trunc`, `svd_full`, `svd_compact`, and `svd_trunc` functions. -Additionally, the default value can also be controlled with a global toggle using [`MatrixAlgebraKit.default_gaugefix`](@ref). +Additionally, the default value can also be controlled with a global toggle using [`MatrixAlgebraKit.default_fixgauge`](@ref). ```@docs; canonical=false MatrixAlgebraKit.gaugefix! -MatrixAlgebraKit.default_gaugefix +MatrixAlgebraKit.default_fixgauge ``` diff --git a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl index e09fccfa..43154b75 100644 --- a/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl +++ b/ext/MatrixAlgebraKitGenericLinearAlgebraExt.jl @@ -1,7 +1,7 @@ module MatrixAlgebraKitGenericLinearAlgebraExt using MatrixAlgebraKit -using MatrixAlgebraKit: sign_safe, check_input, diagview, gaugefix!, default_gaugefix +using MatrixAlgebraKit: sign_safe, check_input, diagview, gaugefix!, default_fixgauge using GenericLinearAlgebra: svd!, svdvals!, eigen!, eigvals!, Hermitian, qr! using LinearAlgebra: I, Diagonal, lmul! @@ -17,7 +17,7 @@ function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix, USVᴴ, alg::GLA_QRIte F = svd!(A) U, S, Vᴴ = F.U, Diagonal(F.S), F.Vt - do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool + do_gauge_fix = get(alg.kwargs, :fixgauge, default_fixgauge())::Bool do_gauge_fix && gaugefix!(svd_compact!, U, Vᴴ) return U, S, Vᴴ @@ -29,7 +29,7 @@ function MatrixAlgebraKit.svd_full!(A::AbstractMatrix, USVᴴ, alg::GLA_QRIterat S = MatrixAlgebraKit.zero!(similar(F.S, (size(U, 2), size(Vᴴ, 1)))) diagview(S) .= F.S - do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool + do_gauge_fix = get(alg.kwargs, :fixgauge, default_fixgauge())::Bool do_gauge_fix && gaugefix!(svd_full!, U, Vᴴ) return U, S, Vᴴ diff --git a/src/common/defaults.jl b/src/common/defaults.jl index 9f48c941..dad16376 100644 --- a/src/common/defaults.jl +++ b/src/common/defaults.jl @@ -43,18 +43,18 @@ Default tolerance for deciding to warn if the provided `A` is not hermitian. default_hermitian_tol(A) = eps(norm(A, Inf))^(3 / 4) -const DEFAULT_GAUGEFIX = Ref(true) +const DEFAULT_FIXGAUGE = Ref(true) @doc """ - default_gaugefix() -> current_value - default_gaugefix(new_value::Bool) -> previous_value + default_fixgauge() -> current_value + default_fixgauge(new_value::Bool) -> previous_value Global toggle for enabling or disabling the default behavior of gauge fixing the output of the eigen- and singular value decompositions. -""" default_gaugefix +""" default_fixgauge -default_gaugefix() = DEFAULT_GAUGEFIX[] -function default_gaugefix(new_value::Bool) - previous_value = DEFAULT_GAUGEFIX[] - DEFAULT_GAUGEFIX[] = new_value +default_fixgauge() = DEFAULT_FIXGAUGE[] +function default_fixgauge(new_value::Bool) + previous_value = DEFAULT_FIXGAUGE[] + DEFAULT_FIXGAUGE[] = new_value return previous_value end diff --git a/src/implementations/eig.jl b/src/implementations/eig.jl index bdb6981a..9b14167c 100644 --- a/src/implementations/eig.jl +++ b/src/implementations/eig.jl @@ -82,8 +82,8 @@ function eig_full!(A::AbstractMatrix, DV, alg::LAPACK_EigAlgorithm) check_input(eig_full!, A, DV, alg) D, V = DV - do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + do_gauge_fix = get(alg.kwargs, :fixgauge, default_fixgauge())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa LAPACK_Simple isempty(alg_kwargs) || @@ -102,7 +102,7 @@ function eig_vals!(A::AbstractMatrix, D, alg::LAPACK_EigAlgorithm) check_input(eig_vals!, A, D, alg) V = similar(A, complex(eltype(A)), (size(A, 1), 0)) - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa LAPACK_Simple isempty(alg_kwargs) || @@ -145,8 +145,8 @@ function eig_full!(A::AbstractMatrix, DV, alg::GPU_EigAlgorithm) check_input(eig_full!, A, DV, alg) D, V = DV - do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + do_gauge_fix = get(alg.kwargs, :fixgauge, default_fixgauge())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa GPU_Simple isempty(alg_kwargs) || @warn "invalid keyword arguments for GPU_Simple" @@ -162,7 +162,7 @@ function eig_vals!(A::AbstractMatrix, D, alg::GPU_EigAlgorithm) check_input(eig_vals!, A, D, alg) V = similar(A, complex(eltype(A)), (size(A, 1), 0)) - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa GPU_Simple isempty(alg_kwargs) || @warn "invalid keyword arguments for GPU_Simple" diff --git a/src/implementations/eigh.jl b/src/implementations/eigh.jl index d4a67dd6..a45300dc 100644 --- a/src/implementations/eigh.jl +++ b/src/implementations/eigh.jl @@ -92,8 +92,8 @@ function eigh_full!(A::AbstractMatrix, DV, alg::LAPACK_EighAlgorithm) D, V = DV Dd = D.diag - do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + do_gauge_fix = get(alg.kwargs, :fixgauge, default_fixgauge())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa LAPACK_MultipleRelativelyRobustRepresentations YALAPACK.heevr!(A, Dd, V; alg_kwargs...) @@ -114,7 +114,7 @@ function eigh_vals!(A::AbstractMatrix, D, alg::LAPACK_EighAlgorithm) check_input(eigh_vals!, A, D, alg) V = similar(A, (size(A, 1), 0)) - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa LAPACK_MultipleRelativelyRobustRepresentations YALAPACK.heevr!(A, D, V; alg_kwargs...) @@ -168,8 +168,8 @@ function eigh_full!(A::AbstractMatrix, DV, alg::GPU_EighAlgorithm) D, V = DV Dd = D.diag - do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + do_gauge_fix = get(alg.kwargs, :fixgauge, default_fixgauge())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa GPU_Jacobi _gpu_heevj!(A, Dd, V; alg_kwargs...) @@ -192,7 +192,7 @@ function eigh_vals!(A::AbstractMatrix, D, alg::GPU_EighAlgorithm) check_input(eigh_vals!, A, D, alg) V = similar(A, (size(A, 1), 0)) - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa GPU_Jacobi _gpu_heevj!(A, D, V; alg_kwargs...) diff --git a/src/implementations/gen_eig.jl b/src/implementations/gen_eig.jl index 94dfe47e..f4a6a7f5 100644 --- a/src/implementations/gen_eig.jl +++ b/src/implementations/gen_eig.jl @@ -58,8 +58,8 @@ function gen_eig_full!(A::AbstractMatrix, B::AbstractMatrix, WV, alg::LAPACK_Eig check_input(gen_eig_full!, A, B, WV, alg) W, V = WV - do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + do_gauge_fix = get(alg.kwargs, :fixgauge, default_fixgauge())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa LAPACK_Simple isempty(alg_kwargs) || @@ -78,7 +78,7 @@ function gen_eig_vals!(A::AbstractMatrix, B::AbstractMatrix, W, alg::LAPACK_EigA check_input(gen_eig_vals!, A, B, W, alg) V = similar(A, complex(eltype(A)), (size(A, 1), 0)) - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa LAPACK_Simple isempty(alg_kwargs) || diff --git a/src/implementations/svd.jl b/src/implementations/svd.jl index 2066b2cd..c8b27f3f 100644 --- a/src/implementations/svd.jl +++ b/src/implementations/svd.jl @@ -120,8 +120,8 @@ function svd_full!(A::AbstractMatrix, USVᴴ, alg::LAPACK_SVDAlgorithm) return USVᴴ end - do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + do_gauge_fix = get(alg.kwargs, :fixgauge, default_fixgauge())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa LAPACK_QRIteration isempty(alg_kwargs) || @@ -153,8 +153,8 @@ function svd_compact!(A::AbstractMatrix, USVᴴ, alg::LAPACK_SVDAlgorithm) check_input(svd_compact!, A, USVᴴ, alg) U, S, Vᴴ = USVᴴ - do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + do_gauge_fix = get(alg.kwargs, :fixgauge, default_fixgauge())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa LAPACK_QRIteration isempty(alg_kwargs) || @@ -183,7 +183,7 @@ function svd_vals!(A::AbstractMatrix, S, alg::LAPACK_SVDAlgorithm) check_input(svd_vals!, A, S, alg) U, Vᴴ = similar(A, (0, 0)), similar(A, (0, 0)) - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa LAPACK_QRIteration isempty(alg_kwargs) || @@ -336,8 +336,8 @@ function svd_full!(A::AbstractMatrix, USVᴴ, alg::GPU_SVDAlgorithm) return USVᴴ end - do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + do_gauge_fix = get(alg.kwargs, :fixgauge, default_fixgauge())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa GPU_QRIteration isempty(alg_kwargs) || @warn "invalid keyword arguments for GPU_QRIteration" @@ -368,7 +368,7 @@ function svd_trunc!(A::AbstractMatrix, USVᴴ, alg::TruncatedAlgorithm{<:GPU_Ran # normal `truncation_error!` does not work here since `S` is not the full singular value spectrum ϵ = sqrt(norm(A)^2 - norm(diagview(Str))^2) # is there a more accurate way to do this? - do_gauge_fix = get(alg.alg.kwargs, :gaugefix, default_gaugefix())::Bool + do_gauge_fix = get(alg.alg.kwargs, :fixgauge, default_fixgauge())::Bool do_gauge_fix && gaugefix!(svd_trunc!, Utr, Vᴴtr) return Utr, Str, Vᴴtr, ϵ @@ -378,8 +378,8 @@ function svd_compact!(A::AbstractMatrix, USVᴴ, alg::GPU_SVDAlgorithm) check_input(svd_compact!, A, USVᴴ, alg) U, S, Vᴴ = USVᴴ - do_gauge_fix = get(alg.kwargs, :gaugefix, default_gaugefix())::Bool - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + do_gauge_fix = get(alg.kwargs, :fixgauge, default_fixgauge())::Bool + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa GPU_QRIteration isempty(alg_kwargs) || @warn "invalid keyword arguments for GPU_QRIteration" @@ -403,7 +403,7 @@ function svd_vals!(A::AbstractMatrix, S, alg::GPU_SVDAlgorithm) check_input(svd_vals!, A, S, alg) U, Vᴴ = similar(A, (0, 0)), similar(A, (0, 0)) - alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:gaugefix,)}) + alg_kwargs = Base.structdiff(alg.kwargs, NamedTuple{(:fixgauge,)}) if alg isa GPU_QRIteration isempty(alg_kwargs) || @warn "invalid keyword arguments for GPU_QRIteration" diff --git a/src/interface/decompositions.jl b/src/interface/decompositions.jl index df5d55e0..0fc7403d 100644 --- a/src/interface/decompositions.jl +++ b/src/interface/decompositions.jl @@ -49,21 +49,21 @@ of `R` are non-negative. # General Eigenvalue Decomposition # ------------------------------- """ - LAPACK_Simple(; gaugefix::Bool = true) + LAPACK_Simple(; fixgauge::Bool = true) Algorithm type to denote the simple LAPACK driver for computing the Schur or non-Hermitian eigenvalue decomposition of a matrix. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigenvectors, +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the eigenvectors, see also [`gaugefix!`](@ref). """ @algdef LAPACK_Simple """ - LAPACK_Expert(; gaugefix::Bool = true) + LAPACK_Expert(; fixgauge::Bool = true) Algorithm type to denote the expert LAPACK driver for computing the Schur or non-Hermitian eigenvalue decomposition of a matrix. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigenvectors, +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the eigenvectors, see also [`gaugefix!`](@ref). """ @algdef LAPACK_Expert @@ -81,44 +81,44 @@ eigenvalue decomposition of a non-Hermitian matrix. # Hermitian Eigenvalue Decomposition # ---------------------------------- """ - LAPACK_QRIteration(; gaugefix::Bool = true) + LAPACK_QRIteration(; fixgauge::Bool = true) Algorithm type to denote the LAPACK driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the QR Iteration algorithm. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the eigen or singular vectors, see also [`gaugefix!`](@ref). """ @algdef LAPACK_QRIteration """ - LAPACK_Bisection(; gaugefix::Bool = true) + LAPACK_Bisection(; fixgauge::Bool = true) Algorithm type to denote the LAPACK driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the Bisection algorithm. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the eigen or singular vectors, see also [`gaugefix!`](@ref). """ @algdef LAPACK_Bisection """ - LAPACK_DivideAndConquer(; gaugefix::Bool = true) + LAPACK_DivideAndConquer(; fixgauge::Bool = true) Algorithm type to denote the LAPACK driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the Divide and Conquer algorithm. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the eigen or singular vectors, see also [`gaugefix!`](@ref). """ @algdef LAPACK_DivideAndConquer """ - LAPACK_MultipleRelativelyRobustRepresentations(; gaugefix::Bool = true) + LAPACK_MultipleRelativelyRobustRepresentations(; fixgauge::Bool = true) Algorithm type to denote the LAPACK driver for computing the eigenvalue decomposition of a Hermitian matrix using the Multiple Relatively Robust Representations algorithm. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigenvectors, +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the eigenvectors, see also [`gaugefix!`](@ref). """ @algdef LAPACK_MultipleRelativelyRobustRepresentations @@ -131,12 +131,12 @@ const LAPACK_EighAlgorithm = Union{ } """ - GLA_QRIteration(; gaugefix::Bool = true) + GLA_QRIteration(; fixgauge::Bool = true) Algorithm type to denote the GenericLinearAlgebra.jl implementation for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the eigen or singular vectors, see also [`gaugefix!`](@ref). """ @algdef GLA_QRIteration @@ -144,11 +144,11 @@ singular vectors, see also [`gaugefix!`](@ref). # Singular Value Decomposition # ---------------------------- """ - LAPACK_Jacobi(; gaugefix::Bool = true) + LAPACK_Jacobi(; fixgauge::Bool = true) Algorithm type to denote the LAPACK driver for computing the singular value decomposition of a general matrix using the Jacobi algorithm. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the singular vectors, +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the singular vectors, see also [`gaugefix!`](@ref). """ @algdef LAPACK_Jacobi @@ -221,33 +221,33 @@ the diagonal elements of `R` are non-negative. @algdef CUSOLVER_HouseholderQR """ - CUSOLVER_QRIteration(; gaugefix::Bool = true) + CUSOLVER_QRIteration(; fixgauge::Bool = true) Algorithm type to denote the CUSOLVER driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the QR Iteration algorithm. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the eigen or singular vectors, see also [`gaugefix!`](@ref). """ @algdef CUSOLVER_QRIteration """ - CUSOLVER_SVDPolar(; gaugefix::Bool = true) + CUSOLVER_SVDPolar(; fixgauge::Bool = true) Algorithm type to denote the CUSOLVER driver for computing the singular value decomposition of a general matrix by using Halley's iterative algorithm to compute the polar decompositon, followed by the hermitian eigenvalue decomposition of the positive definite factor. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the singular +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the singular vectors, see also [`gaugefix!`](@ref). """ @algdef CUSOLVER_SVDPolar """ - CUSOLVER_Jacobi(; gaugefix::Bool = true) + CUSOLVER_Jacobi(; fixgauge::Bool = true) Algorithm type to denote the CUSOLVER driver for computing the singular value decomposition of a general matrix using the Jacobi algorithm. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the singular +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the singular vectors, see also [`gaugefix!`](@ref). """ @algdef CUSOLVER_Jacobi @@ -270,11 +270,11 @@ for more information. does_truncate(::TruncatedAlgorithm{<:CUSOLVER_Randomized}) = true """ - CUSOLVER_Simple(; gaugefix::Bool = true) + CUSOLVER_Simple(; fixgauge::Bool = true) Algorithm type to denote the simple CUSOLVER driver for computing the non-Hermitian eigenvalue decomposition of a matrix. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigenvectors, +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the eigenvectors, see also [`gaugefix!`](@ref). """ @algdef CUSOLVER_Simple @@ -282,12 +282,12 @@ see also [`gaugefix!`](@ref). const CUSOLVER_EigAlgorithm = Union{CUSOLVER_Simple} """ - CUSOLVER_DivideAndConquer(; gaugefix::Bool = true) + CUSOLVER_DivideAndConquer(; fixgauge::Bool = true) Algorithm type to denote the CUSOLVER driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the Divide and Conquer algorithm. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the eigen or singular vectors, see also [`gaugefix!`](@ref). """ @algdef CUSOLVER_DivideAndConquer @@ -309,44 +309,44 @@ the diagonal elements of `R` are non-negative. @algdef ROCSOLVER_HouseholderQR """ - ROCSOLVER_QRIteration(; gaugefix::Bool = true) + ROCSOLVER_QRIteration(; fixgauge::Bool = true) Algorithm type to denote the ROCSOLVER driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the QR Iteration algorithm. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the eigen or singular vectors, see also [`gaugefix!`](@ref). """ @algdef ROCSOLVER_QRIteration """ - ROCSOLVER_Jacobi(; gaugefix::Bool = true) + ROCSOLVER_Jacobi(; fixgauge::Bool = true) Algorithm type to denote the ROCSOLVER driver for computing the singular value decomposition of a general matrix using the Jacobi algorithm. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the singular +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the singular vectors, see also [`gaugefix!`](@ref). """ @algdef ROCSOLVER_Jacobi """ - ROCSOLVER_Bisection(; gaugefix::Bool = true) + ROCSOLVER_Bisection(; fixgauge::Bool = true) Algorithm type to denote the ROCSOLVER driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the Bisection algorithm. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the eigen or singular vectors, see also [`gaugefix!`](@ref). """ @algdef ROCSOLVER_Bisection """ - ROCSOLVER_DivideAndConquer(; gaugefix::Bool = true) + ROCSOLVER_DivideAndConquer(; fixgauge::Bool = true) Algorithm type to denote the ROCSOLVER driver for computing the eigenvalue decomposition of a Hermitian matrix, or the singular value decomposition of a general matrix using the Divide and Conquer algorithm. -The `gaugefix` keyword can be used to toggle whether or not to fix the gauge of the eigen or +The `fixgauge` keyword can be used to toggle whether or not to fix the gauge of the eigen or singular vectors, see also [`gaugefix!`](@ref). """ @algdef ROCSOLVER_DivideAndConquer diff --git a/src/pullbacks/eig.jl b/src/pullbacks/eig.jl index 3115b3d5..ff0de512 100644 --- a/src/pullbacks/eig.jl +++ b/src/pullbacks/eig.jl @@ -43,7 +43,7 @@ function eig_pullback!( mask = abs.(transpose(D) .- D) .< degeneracy_atol Δgauge = norm(view(VᴴΔV, mask), Inf) - Δgauge < gauge_atol || + Δgauge ≤ gauge_atol || @warn "`eig` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" VᴴΔV .*= conj.(inv_safe.(transpose(D) .- D, degeneracy_atol)) @@ -121,7 +121,7 @@ function eig_trunc_pullback!( VᴴΔV = V' * ΔV mask = abs.(transpose(D) .- D) .< degeneracy_atol Δgauge = norm(view(VᴴΔV, mask), Inf) - Δgauge < gauge_atol || + Δgauge ≤ gauge_atol || @warn "`eig` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" ΔVperp = ΔV - V * inv(G) * VᴴΔV diff --git a/src/pullbacks/eigh.jl b/src/pullbacks/eigh.jl index e7c93c75..d303faa7 100644 --- a/src/pullbacks/eigh.jl +++ b/src/pullbacks/eigh.jl @@ -44,7 +44,7 @@ function eigh_pullback!( mask = abs.(D' .- D) .< degeneracy_atol Δgauge = norm(view(aVᴴΔV, mask)) - Δgauge < gauge_atol || + Δgauge ≤ gauge_atol || @warn "`eigh` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" aVᴴΔV .*= inv_safe.(D' .- D, degeneracy_atol) @@ -112,7 +112,7 @@ function eigh_trunc_pullback!( mask = abs.(D' .- D) .< degeneracy_atol Δgauge = norm(view(aVᴴΔV, mask)) - Δgauge < gauge_atol || + Δgauge ≤ gauge_atol || @warn "`eigh` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" aVᴴΔV .*= inv_safe.(D' .- D, degeneracy_atol) diff --git a/src/pullbacks/lq.jl b/src/pullbacks/lq.jl index d2cfa290..b30fe198 100644 --- a/src/pullbacks/lq.jl +++ b/src/pullbacks/lq.jl @@ -50,7 +50,7 @@ function lq_pullback!( ΔL22 = view(ΔL, (p + 1):m, (p + 1):minmn) Δgauge = max(Δgauge, norm(ΔL22, Inf)) end - Δgauge < gauge_atol || + Δgauge ≤ gauge_atol || @warn "`lq` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" end @@ -70,7 +70,7 @@ function lq_pullback!( # Q2' * ΔQ2 as a gauge dependent quantity. ΔQ2Q1ᴴ = ΔQ2 * Q1' Δgauge = norm(mul!(copy(ΔQ2), ΔQ2Q1ᴴ, Q1, -1, 1), Inf) - Δgauge < gauge_atol || + Δgauge ≤ gauge_atol || @warn "`lq` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" ΔQ̃ = mul!(ΔQ̃, ΔQ2Q1ᴴ', Q2, -1, 1) end @@ -120,7 +120,7 @@ function lq_null_pullback!( if !iszerotangent(ΔNᴴ) && size(Nᴴ, 1) > 0 aNᴴΔN = project_antihermitian!(Nᴴ * ΔNᴴ') Δgauge = norm(aNᴴΔN) - Δgauge < gauge_atol || + Δgauge ≤ gauge_atol || @warn "`lq_null` cotangent sensitive to gauge choice: (|Δgauge| = $Δgauge)" L, Q = lq_compact(A; positive = true) # should we be able to provide algorithm here? X = ldiv!(LowerTriangular(L)', Q * ΔNᴴ') diff --git a/src/pullbacks/qr.jl b/src/pullbacks/qr.jl index 10882bdf..888029be 100644 --- a/src/pullbacks/qr.jl +++ b/src/pullbacks/qr.jl @@ -51,7 +51,7 @@ function qr_pullback!( ΔR22 = view(ΔR, (p + 1):minmn, (p + 1):n) Δgauge = max(Δgauge, norm(ΔR22, Inf)) end - Δgauge < gauge_atol || + Δgauge ≤ gauge_atol || @warn "`qr` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" end @@ -70,7 +70,7 @@ function qr_pullback!( # Q2' * ΔQ2 as a gauge dependent quantity. Q1dΔQ2 = Q1' * ΔQ2 Δgauge = norm(mul!(copy(ΔQ2), Q1, Q1dΔQ2, -1, 1), Inf) - Δgauge < gauge_atol || + Δgauge ≤ gauge_atol || @warn "`qr` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" ΔQ̃ = mul!(ΔQ̃, Q2, Q1dΔQ2', -1, 1) end @@ -120,7 +120,7 @@ function qr_null_pullback!( if !iszerotangent(ΔN) && size(N, 2) > 0 aNᴴΔN = project_antihermitian!(N' * ΔN) Δgauge = norm(aNᴴΔN) - Δgauge < gauge_atol || + Δgauge ≤ gauge_atol || @warn "`qr_null` cotangent sensitive to gauge choice: (|Δgauge| = $Δgauge)" Q, R = qr_compact(A; positive = true) diff --git a/src/pullbacks/svd.jl b/src/pullbacks/svd.jl index a85b7165..1667be65 100644 --- a/src/pullbacks/svd.jl +++ b/src/pullbacks/svd.jl @@ -72,7 +72,7 @@ function svd_pullback!( # check whether cotangents arise from gauge-invariance objective function mask = abs.(Sr' .- Sr) .< degeneracy_atol Δgauge = norm(view(aUΔU, mask) + view(aVΔV, mask), Inf) - Δgauge < gauge_atol || + Δgauge ≤ gauge_atol || @warn "`svd` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" UdΔAV = (aUΔU .+ aVΔV) .* inv_safe.(Sr' .- Sr, degeneracy_atol) .+ @@ -160,7 +160,7 @@ function svd_trunc_pullback!( # check whether cotangents arise from gauge-invariance objective function mask = abs.(S' .- S) .< degeneracy_atol Δgauge = norm(view(aUΔU, mask) + view(aVΔV, mask), Inf) - Δgauge < gauge_atol || + Δgauge ≤ gauge_atol || @warn "`svd` cotangents sensitive to gauge choice: (|Δgauge| = $Δgauge)" UdΔAV = (aUΔU .+ aVΔV) .* inv_safe.(S' .- S, degeneracy_atol) .+ diff --git a/test/mooncake.jl b/test/mooncake.jl index 554b35d9..7aed68cf 100644 --- a/test/mooncake.jl +++ b/test/mooncake.jl @@ -504,6 +504,20 @@ end end end +left_orth_qr(X) = left_orth(X; alg = :qr) +left_orth_polar(X) = left_orth(X; alg = :polar) +left_null_qr(X) = left_null(X; alg = :qr) +right_orth_lq(X) = right_orth(X; alg = :lq) +right_orth_polar(X) = right_orth(X; alg = :polar) +right_null_lq(X) = right_null(X; alg = :lq) + +MatrixAlgebraKit.copy_input(::typeof(left_orth_qr), A) = MatrixAlgebraKit.copy_input(left_orth, A) +MatrixAlgebraKit.copy_input(::typeof(left_orth_polar), A) = MatrixAlgebraKit.copy_input(left_orth, A) +MatrixAlgebraKit.copy_input(::typeof(left_null_qr), A) = MatrixAlgebraKit.copy_input(left_null, A) +MatrixAlgebraKit.copy_input(::typeof(right_orth_lq), A) = MatrixAlgebraKit.copy_input(right_orth, A) +MatrixAlgebraKit.copy_input(::typeof(right_orth_polar), A) = MatrixAlgebraKit.copy_input(right_orth, A) +MatrixAlgebraKit.copy_input(::typeof(right_null_lq), A) = MatrixAlgebraKit.copy_input(right_null, A) + @timedtestset "Orth and null with eltype $T" for T in ETs rng = StableRNG(12345) m = 19 @@ -517,11 +531,6 @@ end Mooncake.TestUtils.test_rule(rng, right_orth, A; mode = Mooncake.ReverseMode, atol = atol, rtol = rtol, is_primitive = false) test_pullbacks_match(rng, right_orth!, right_orth, A, CVᴴ, (randn(rng, T, size(CVᴴ[1])...), randn(rng, T, size(CVᴴ[2])...))) - left_orth_qr(X) = left_orth(X; alg = :qr) - left_orth_polar(X) = left_orth(X; alg = :polar) - MatrixAlgebraKit.copy_input(left_orth_qr, A) = MatrixAlgebraKit.copy_input(left_orth, A) - MatrixAlgebraKit.copy_input(left_orth_polar, A) = MatrixAlgebraKit.copy_input(left_orth, A) - Mooncake.TestUtils.test_rule(rng, left_orth_qr, A; mode = Mooncake.ReverseMode, atol = atol, rtol = rtol, is_primitive = false) test_pullbacks_match(rng, ((X, VC) -> left_orth!(X, VC; alg = :qr)), left_orth_qr, A, VC, (randn(rng, T, size(VC[1])...), randn(rng, T, size(VC[2])...))) if m >= n @@ -529,18 +538,12 @@ end test_pullbacks_match(rng, ((X, VC) -> left_orth!(X, VC; alg = :polar)), left_orth_polar, A, VC, (randn(rng, T, size(VC[1])...), randn(rng, T, size(VC[2])...))) end - left_null_qr(X) = left_null(X; alg = :qr) - MatrixAlgebraKit.copy_input(left_null_qr, A) = MatrixAlgebraKit.copy_input(left_null, A) N = left_orth(A; alg = :qr)[1] * randn(rng, T, min(m, n), m - min(m, n)) ΔN = left_orth(A; alg = :qr)[1] * randn(rng, T, min(m, n), m - min(m, n)) dN = make_mooncake_tangent(ΔN) Mooncake.TestUtils.test_rule(rng, left_null_qr, A; mode = Mooncake.ReverseMode, atol = atol, rtol = rtol, is_primitive = false, output_tangent = dN) test_pullbacks_match(rng, ((X, N) -> left_null!(X, N; alg = :qr)), left_null_qr, A, N, ΔN) - right_orth_lq(X) = right_orth(X; alg = :lq) - right_orth_polar(X) = right_orth(X; alg = :polar) - MatrixAlgebraKit.copy_input(right_orth_lq, A) = MatrixAlgebraKit.copy_input(right_orth, A) - MatrixAlgebraKit.copy_input(right_orth_polar, A) = MatrixAlgebraKit.copy_input(right_orth, A) Mooncake.TestUtils.test_rule(rng, right_orth_lq, A; mode = Mooncake.ReverseMode, atol = atol, rtol = rtol, is_primitive = false) test_pullbacks_match(rng, ((X, CVᴴ) -> right_orth!(X, CVᴴ; alg = :lq)), right_orth_lq, A, CVᴴ, (randn(rng, T, size(CVᴴ[1])...), randn(rng, T, size(CVᴴ[2])...))) @@ -549,8 +552,6 @@ end test_pullbacks_match(rng, ((X, CVᴴ) -> right_orth!(X, CVᴴ; alg = :polar)), right_orth_polar, A, CVᴴ, (randn(rng, T, size(CVᴴ[1])...), randn(rng, T, size(CVᴴ[2])...))) end - right_null_lq(X) = right_null(X; alg = :lq) - MatrixAlgebraKit.copy_input(right_null_lq, A) = MatrixAlgebraKit.copy_input(right_null, A) Nᴴ = randn(rng, T, n - min(m, n), min(m, n)) * right_orth(A; alg = :lq)[2] ΔNᴴ = randn(rng, T, n - min(m, n), min(m, n)) * right_orth(A; alg = :lq)[2] dNᴴ = make_mooncake_tangent(ΔNᴴ)