some final fixes

Author: Jutho

docs/src/user_interface/truncations.md

@@ -23,3 +23,17 @@ maxdim = 10
 atol = 1e-6
 combined_trunc = truncrank(maxdim) & trunctol(; atol)
 ```
+
+## Truncation Error
+
+When using truncated decompositions such as [`svd_trunc`](@ref), [`eig_trunc`](@ref), or [`eigh_trunc`](@ref), an additional truncation error value is returned.
+This error is defined as the 2-norm of the discarded  singular values or eigenvalues, providing a measure of the approximation quality.
+For `svd_trunc` and `eigh_trunc`, this corresponds to the 2-norm difference between the original and the truncated matrix.
+For the case of `eig_trunc`, this interpretation does not hold because the norm of the non-unitary matrix of eigenvectors and its inverse also influence the approximation quality.
+
+
+For example:
+```julia
+U, S, Vᴴ, ϵ = svd_trunc(A; trunc=truncrank(10))
+# ϵ is the 2-norm of the discarded singular values
+```

ext/MatrixAlgebraKitChainRulesCoreExt.jl

@@ -113,7 +113,7 @@ for eig in (:eig, :eigh)
             Ac = copy_input($eig_f, A)
             DV = $(eig_f!)(Ac, DV, alg.alg)
             DV′, ind = MatrixAlgebraKit.truncate($eig_t!, DV, alg.trunc)
-            ϵ = compute_truncerr!(diagview(copy(DV[1])), ind)
+            ϵ = compute_truncerr!(copy(diagview(DV[1])), ind)
             return (DV′..., ϵ), $(_make_eig_t_pb)(A, DV, ind)
         end
         function $(_make_eig_t_pb)(A, DV, ind)
@@ -157,7 +157,7 @@ function ChainRulesCore.rrule(::typeof(svd_trunc!), A, USVᴴ, alg::TruncatedAlg
     Ac = copy_input(svd_compact, A)
     USVᴴ = svd_compact!(Ac, USVᴴ, alg.alg)
     USVᴴ′, ind = MatrixAlgebraKit.truncate(svd_trunc!, USVᴴ, alg.trunc)
-    ϵ = compute_truncerr!(diagview(copy(USVᴴ[2])), ind)
+    ϵ = compute_truncerr!(copy(diagview(USVᴴ[2])), ind)
     return (USVᴴ′..., ϵ), _make_svd_trunc_pullback(A, USVᴴ, ind)
 end
 function _make_svd_trunc_pullback(A, USVᴴ, ind)

src/implementations/svd.jl

@@ -382,7 +382,12 @@ function svd_trunc!(A::AbstractMatrix, USVᴴ, alg::TruncatedAlgorithm{<:GPU_Ran
     _gpu_Xgesvdr!(A, S.diag, U, Vᴴ; alg.alg.kwargs...)
     # TODO: make this controllable using a `gaugefix` keyword argument
     gaugefix!(svd_trunc!, U, S, Vᴴ, size(A)...)
-    return first(truncate(svd_trunc!, USVᴴ, alg.trunc))
+    # TODO: make sure that truncation is based on maxrank, otherwise this might be wrong
+    USVᴴtrunc, ind = truncate(svd_trunc!, (U, S, Vᴴ), alg.trunc)
+    Strunc = diagview(USVᴴtrunc[2])
+    # normal `compute_truncerr!` does not work here since `S` is not the full singular value spectrum
+    ϵ = sqrt(norm(A)^2 - norm(Strunc)^2) # is there a more accurate way to do this?
+    return USVᴴtrunc..., ϵ
 end
 
 function MatrixAlgebraKit.svd_compact!(A::AbstractMatrix, USVᴴ, alg::GPU_SVDAlgorithm)

test/cuda/svd.jl

@@ -106,6 +106,7 @@ end
             U1, S1, V1ᴴ, ϵ1 = @constinferred svd_trunc(A; alg, trunc = truncrank(r))
             @test length(S1.diag) == r
             @test opnorm(A - U1 * S1 * V1ᴴ) ≈ S₀[r + 1]
+            @test norm(A - U1 * S1 * V1ᴴ) ≈ ϵ1
 
             if !(alg isa CUSOLVER_Randomized)
                 s = 1 + sqrt(eps(real(T)))