From a47d84f15fdaf783dd6f18f608fcd14f0675e9eb Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Fri, 12 Jun 2026 14:59:56 -0500 Subject: [PATCH] Apply runic formatting --- src/GsvdInitialization.jl | 102 ++++++++++++++++++----------------- test/runtests.jl | 110 +++++++++++++++++++------------------- 2 files changed, 108 insertions(+), 104 deletions(-) diff --git a/src/GsvdInitialization.jl b/src/GsvdInitialization.jl index 77b2cc7..e7954fb 100644 --- a/src/GsvdInitialization.jl +++ b/src/GsvdInitialization.jl @@ -19,7 +19,7 @@ using Kronecker: kronecker using SparseArrays: sparse export gsvdnmf, - gsvdrecover + gsvdrecover @static if VERSION >= v"1.11" eval(Meta.parse("public truncating, joint_nnls")) @@ -83,12 +83,14 @@ julia> size(result.W), size(result.H) ((4, 3), (3, 5)) ``` """ -function gsvdnmf(strategy, X::AbstractMatrix, W::AbstractMatrix, H::AbstractMatrix, f; - n2 = size(first(f), 2), - tol_final = 1e-4, - alg = :cd, - truncmult = 1e-5, - kwargs...) +function gsvdnmf( + strategy, X::AbstractMatrix, W::AbstractMatrix, H::AbstractMatrix, f; + n2 = size(first(f), 2), + tol_final = 1.0e-4, + alg = :cd, + truncmult = 1.0e-5, + kwargs... + ) Base.require_one_based_indexing(X, W, H) n1 = size(W, 2) kadd = n2 - n1 @@ -99,7 +101,7 @@ function gsvdnmf(strategy, X::AbstractMatrix, W::AbstractMatrix, H::AbstractMatr if alg == :multmse W_recover, H_recover = max.(W_recover, truncmult), max.(H_recover, truncmult) end - result_recover = nnmf(X, n2; kwargs..., alg, init=:custom, tol=tol_final, W0=copy(W_recover), H0=copy(H_recover)) + result_recover = nnmf(X, n2; kwargs..., alg, init = :custom, tol = tol_final, W0 = copy(W_recover), H0 = copy(H_recover)) return result_recover, Λ end gsvdnmf(X::AbstractMatrix, W::AbstractMatrix, H::AbstractMatrix, f; kwargs...) = @@ -156,19 +158,19 @@ julia> sum(abs2, X - result.W*result.H) < 1e-6 * sum(abs2, X) # near-exact rank true ``` """ -function gsvdnmf(strategy, X::AbstractMatrix, ncomponents::Pair{<:Integer,<:Integer}; tol_final=1e-4, tol_intermediate=tol_final, kwargs...) +function gsvdnmf(strategy, X::AbstractMatrix, ncomponents::Pair{<:Integer, <:Integer}; tol_final = 1.0e-4, tol_intermediate = tol_final, kwargs...) Base.require_one_based_indexing(X) n1, n2 = ncomponents f = tsvd(X, n2) W0, H0 = nndsvd(X, n1; initdata = (U = f[1], S = f[2], V = f[3])) - result_initial_nmf = nnmf(X, n1; kwargs..., init=:custom, tol=tol_intermediate, W0=copy(W0), H0=copy(H0)) + result_initial_nmf = nnmf(X, n1; kwargs..., init = :custom, tol = tol_intermediate, W0 = copy(W0), H0 = copy(H0)) W_initial_nmf, H_initial_nmf = result_initial_nmf.W, result_initial_nmf.H return gsvdnmf(strategy, X, W_initial_nmf, H_initial_nmf, f; kwargs..., n2, tol_final) end -gsvdnmf(X::AbstractMatrix, ncomponents::Pair{<:Integer,<:Integer}; kwargs...) = +gsvdnmf(X::AbstractMatrix, ncomponents::Pair{<:Integer, <:Integer}; kwargs...) = gsvdnmf(truncating, X, ncomponents; kwargs...) gsvdnmf(strategy, X::AbstractMatrix, ncomponents_final::Integer; kwargs...) = - gsvdnmf(strategy, X, ncomponents_final-1 => ncomponents_final; kwargs...) + gsvdnmf(strategy, X, ncomponents_final - 1 => ncomponents_final; kwargs...) gsvdnmf(X::AbstractMatrix, ncomponents_final::Integer; kwargs...) = gsvdnmf(truncating, X, ncomponents_final; kwargs...) @@ -238,7 +240,7 @@ function gsvdrecover(strategy, X, W0::AbstractMatrix, H0::AbstractMatrix, kadd:: # An offset-axes SVD wider than `n` would make the `1:n` slices below # succeed on the wrong columns; reject it before slicing. Base.require_one_based_indexing(U0, S0, V0) - U0, S0, V0 = U0[:,1:n], S0[1:n], V0[:,1:n] + U0, S0, V0 = U0[:, 1:n], S0[1:n], V0[:, 1:n] Hadd, Λ = init_H(U0, S0, V0, W0, H0, kadd) W, H = strategy(X, W0, H0, Hadd) return W, H, Λ @@ -300,12 +302,12 @@ function joint_nnls(X, W0::AbstractMatrix, H0::AbstractMatrix, Hadd::AbstractMat end function init_H(U0::AbstractMatrix, S0::AbstractVector, V0::AbstractMatrix, W0::AbstractMatrix, H0::AbstractMatrix, kadd::Integer) - _, _, Q, D1, D2, R = svd(Matrix(Diagonal(S0)), (U0'*W0)*(H0*V0)); + _, _, Q, D1, D2, R = svd(Matrix(Diagonal(S0)), (U0' * W0) * (H0 * V0)) r0 = size(U0, 2) - k = findfirst(x->x!=0, D2[1,:]) - k = (k === nothing) ? r0 : k-1 + k = findfirst(x -> x != 0, D2[1, :]) + k = (k === nothing) ? r0 : k - 1 kadd >= k || @warn "kadd ($kadd) is less than the rank deficiency of W0*H0 ($k)." - F = (diag(D1[k+1:r0, k+1:r0])./diag(D2[1:r0-k,k+1:r0])).^2 + F = (diag(D1[(k + 1):r0, (k + 1):r0]) ./ diag(D2[1:(r0 - k), (k + 1):r0])) .^ 2 Λ = vcat(fill(Inf, k), F) H_index = sortperm(Λ, rev = true)[1:kadd] # Columns of inv(R*Q') = Q*inv(R) selected by H_index, via a triangular @@ -316,32 +318,32 @@ function init_H(U0::AbstractMatrix, S0::AbstractVector, V0::AbstractMatrix, W0:: E[idx, j] = 1 end Hadd = Q * (UpperTriangular(R) \ E) - Hadd_1 = V0*Hadd + Hadd_1 = V0 * Hadd return Hadd_1', Λ end -function init_W_joint_nnls(X, W0::AbstractMatrix{T}, H0::AbstractMatrix{T}, Hadd::AbstractMatrix{T}) where T +function init_W_joint_nnls(X, W0::AbstractMatrix{T}, H0::AbstractMatrix{T}, Hadd::AbstractMatrix{T}) where {T} m = size(X, 1) kadd = size(Hadd, 1) G = gram_sp_C(W0, H0, Hadd)[1] b = gram_b(X, W0, H0, Hadd) - θ = nonneg_lsq(G, b; alg=:fnnls, gram=true) - Wadd = reshape(θ[1:m*kadd], m, kadd) - α = θ[m*kadd+1:end] + θ = nonneg_lsq(G, b; alg = :fnnls, gram = true) + Wadd = reshape(θ[1:(m * kadd)], m, kadd) + α = θ[(m * kadd + 1):end] return Wadd, α end function gram_sp_C(W0, H0, Hadd) m, r0 = size(W0) k = size(Hadd, 1) - mk = m*k - W0W0, H0H0 = W0'*W0, H0*H0' - P = Hadd*H0' - HH = Hadd*Hadd' - G22 = sparse(W0W0.*H0H0) + mk = m * k + W0W0, H0H0 = W0' * W0, H0 * H0' + P = Hadd * H0' + HH = Hadd * Hadd' + G22 = sparse(W0W0 .* H0H0) G12 = zeros(eltype(W0W0), mk, r0) for j in 1:r0 - G12[:,j] .= vec(W0[:,j] * P[:,j]') + G12[:, j] .= vec(W0[:, j] * P[:, j]') end G12 = sparse(G12) G11 = kronecker(HH, sparse(I, m, m)) @@ -354,49 +356,49 @@ function gram_b(X, W0, H0, Hadd) return b end -function init_W(X, W0::AbstractMatrix{T}, H0::AbstractMatrix{T}, Hadd::AbstractMatrix{T}; α = nothing) where T +function init_W(X, W0::AbstractMatrix{T}, H0::AbstractMatrix{T}, Hadd::AbstractMatrix{T}; α = nothing) where {T} A, b, _, cholHH, H0Hadd, XHaddt = obj_para(X, W0, H0, Hadd) if α === nothing if isposdef(A) - α = nonneg_lsq(A, -b; alg=:fnnls, gram=true) + α = nonneg_lsq(A, -b; alg = :fnnls, gram = true) else # A is not positive definite: the QP min_{α≥0} α'Aα + 2b'α has no # unique bounded minimum, so fnnls is not meaningful. Fall back to # α = 1 (keep existing components at their current scale). - sum(abs2, A) <= 1e-12 || @warn "A is not positive definite." maxlog=1 + sum(abs2, A) <= 1.0e-12 || @warn "A is not positive definite." maxlog = 1 α = ones(T, size(A, 1)) end end - Wadd = (XHaddt - W0*Diagonal(α[:])*H0Hadd) / cholHH + Wadd = (XHaddt - W0 * Diagonal(α[:]) * H0Hadd) / cholHH return Wadd, abs.(α) end -function obj_para(X, W0::AbstractMatrix{T}, H0::AbstractMatrix{T}, Hadd::AbstractMatrix{T}) where T - XHaddt = X*Hadd' - H0Hadd = H0*Hadd' - HH = Hadd*Hadd' - W0W0 = W0'*W0 - H0H0 = H0*H0' +function obj_para(X, W0::AbstractMatrix{T}, H0::AbstractMatrix{T}, Hadd::AbstractMatrix{T}) where {T} + XHaddt = X * Hadd' + H0Hadd = H0 * Hadd' + HH = Hadd * Hadd' + W0W0 = W0' * W0 + H0H0 = H0 * H0' cholHH = cholesky(Symmetric(HH)) - A = W0W0.*(H0H0-H0Hadd*(cholHH \ H0Hadd')) - W0tXH0t = W0'*X*H0' - W0XHaddt = W0'*XHaddt - b = diag(H0Hadd*(cholHH \ W0XHaddt')-W0tXH0t) - C = sum(abs2, X)-tr(cholHH \ (XHaddt'*XHaddt)) + A = W0W0 .* (H0H0 - H0Hadd * (cholHH \ H0Hadd')) + W0tXH0t = W0' * X * H0' + W0XHaddt = W0' * XHaddt + b = diag(H0Hadd * (cholHH \ W0XHaddt') - W0tXH0t) + C = sum(abs2, X) - tr(cholHH \ (XHaddt' * XHaddt)) return Symmetric(A), b, C, cholHH, H0Hadd, XHaddt end -function Wcols_modification(X, W::AbstractMatrix{T}, H::AbstractMatrix{T}) where T - WW, HH = W'*W, H*H' - WtXHt = W'*X*H' +function Wcols_modification(X, W::AbstractMatrix{T}, H::AbstractMatrix{T}) where {T} + WW, HH = W' * W, H * H' + WtXHt = W' * X * H' a = diag(WtXHt) - B = WW.*HH - β = nonneg_lsq(B, a; alg=:fnnls, gram=true) + B = WW .* HH + β = nonneg_lsq(B, a; alg = :fnnls, gram = true) return β[:] end function truncatepos(Y, X, W, H) - ΔX = max.(zero(eltype(X)), X - W*H) + ΔX = max.(zero(eltype(X)), X - W * H) Yout = similar(Y) for j in axes(Y, 2) y = view(Y, :, j) @@ -412,4 +414,4 @@ function truncatepos(Y, X, W, H) end -end \ No newline at end of file +end diff --git a/test/runtests.jl b/test/runtests.jl index 1db5900..e8335cf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -14,13 +14,15 @@ end @testset "ExplicitImports" begin # `nndsvd` is exported by NMF but not declared `public` in NMF.jl. - test_explicit_imports(GsvdInitialization; - ignore = (:nndsvd,), - all_explicit_imports_are_public = VERSION >= v"1.11", - all_qualified_accesses_are_public = VERSION >= v"1.11") + test_explicit_imports( + GsvdInitialization; + ignore = (:nndsvd,), + all_explicit_imports_are_public = VERSION >= v"1.11", + all_qualified_accesses_are_public = VERSION >= v"1.11" + ) end -DocMeta.setdocmeta!(GsvdInitialization, :DocTestSetup, :(using GsvdInitialization); recursive=true) +DocMeta.setdocmeta!(GsvdInitialization, :DocTestSetup, :(using GsvdInitialization); recursive = true) @testset "Doctests" begin doctest(GsvdInitialization; manual = false) end @@ -55,13 +57,13 @@ svdX = load_svd_of_gt() rng = StableRNG(1) W = W_GT H = H_GT - X = W*H - standard_nmf = nnmf(X, 10; alg = :cd, init=:nndsvd, tol=1e-4, maxiter = 10^5, initdata = svdX) - result_gsvd, Λ_gsvd = gsvdnmf(X, 9=>10; alg = :cd, maxiter = 10^5, tol_final=1e-4, tol_intermediate = 1e-4); + X = W * H + standard_nmf = nnmf(X, 10; alg = :cd, init = :nndsvd, tol = 1.0e-4, maxiter = 10^5, initdata = svdX) + result_gsvd, Λ_gsvd = gsvdnmf(X, 9 => 10; alg = :cd, maxiter = 10^5, tol_final = 1.0e-4, tol_intermediate = 1.0e-4) W_gsvd, H_gsvd = result_gsvd.W, result_gsvd.H @test size(W_gsvd, 2) == 10 - @test sum(abs2, X-W_gsvd*H_gsvd)/sum(abs2, X) < 2e-10 - @test sum(abs2, X-standard_nmf.W*standard_nmf.H)/sum(abs2, X) > sum(abs2, X-W_gsvd*H_gsvd)/sum(abs2, X) + @test sum(abs2, X - W_gsvd * H_gsvd) / sum(abs2, X) < 2.0e-10 + @test sum(abs2, X - standard_nmf.W * standard_nmf.H) / sum(abs2, X) > sum(abs2, X - W_gsvd * H_gsvd) / sum(abs2, X) @test length(Λ_gsvd) == 9 # `gsvdnmf(X, n2)` is sugar for `gsvdnmf(X, n2-1 => n2)`. The two calls run @@ -69,12 +71,12 @@ svdX = load_svd_of_gt() # multithreaded BLAS in `nnmf`; `rtol = 1e-6` stays far tighter than the # `1e-4` NMF tol. X = rand(rng, 30, 20) - result_1, _ = gsvdnmf(X, 10; alg=:cd) - result_2, _ = gsvdnmf(X, 9 => 10; alg=:cd) + result_1, _ = gsvdnmf(X, 10; alg = :cd) + result_2, _ = gsvdnmf(X, 9 => 10; alg = :cd) W_gsvd_1, H_gsvd_1 = result_1.W, result_1.H W_gsvd_2, H_gsvd_2 = result_2.W, result_2.H - @test isapprox(W_gsvd_1, W_gsvd_2; rtol = 1e-6) - @test isapprox(H_gsvd_1, H_gsvd_2; rtol = 1e-6) + @test isapprox(W_gsvd_1, W_gsvd_2; rtol = 1.0e-6) + @test isapprox(H_gsvd_1, H_gsvd_2; rtol = 1.0e-6) # n2 == size(W, 2) is a caller bug: there is nothing to augment. Reject it # eagerly rather than silently returning the input factorization. @@ -103,38 +105,38 @@ end @testset "GsvdInitialization" begin rng = StableRNG(2) W, H = rand(rng, 10, 3), rand(rng, 3, 8) - X = W*H + X = W * H U, S, V = svd(X) W0, H0 = copy(W), copy(H) Hadd = rand(rng, 2, 8) Wadd, a = GsvdInitialization.init_W(X, W0, H0, Hadd) @test a ≈ ones(size(W0, 2)) - @test norm(Wadd) <= 1e-8 + @test norm(Wadd) <= 1.0e-8 W0, H0 = zero(W), zero(H) - Hadd = V[:,1:3]' + Hadd = V[:, 1:3]' Wadd, a = GsvdInitialization.init_W(X, W0, H0, Hadd) - @test sum(abs2, Wadd-(U*Diagonal(S))[:,1:3]) <= 1e-12 + @test sum(abs2, Wadd - (U * Diagonal(S))[:, 1:3]) <= 1.0e-12 W0, H0 = rand(rng, 10, 4), rand(rng, 4, 8) Hadd = rand(rng, 2, 8) A, b, C, HH, γ = GsvdInitialization.obj_para(X, W0, H0, Hadd) a = rand(rng, 4) Wadd, a = GsvdInitialization.init_W(X, W0, H0, Hadd, α = a) - E = a'*A*a+2*b'*a+C - @test abs(E-sum(abs2, X-[repeat(a', size(W0, 1)).*W0 Wadd]*[H0;Hadd])) <= 1e-12 + E = a' * A * a + 2 * b' * a + C + @test abs(E - sum(abs2, X - [repeat(a', size(W0, 1)) .* W0 Wadd] * [H0;Hadd])) <= 1.0e-12 β0 = rand(rng, 3) - β = GsvdInitialization.Wcols_modification(X, repeat(β0', size(W, 1)).*W, H) - @test β.*β0 ≈ ones(3) + β = GsvdInitialization.Wcols_modification(X, repeat(β0', size(W, 1)) .* W, H) + @test β .* β0 ≈ ones(3) # When H0 is parallel to Hadd the Schur complement vanishes and A = 0, making # the QP degenerate. init_W must return finite results rather than throwing # SingularException (which fnnls raises on Julia ≥ 1.12 for a zero pivot). W0_deg = rand(rng, Float64, 5, 1) H0_deg = rand(rng, Float64, 1, 8) - X_deg = W0_deg * H0_deg + X_deg = W0_deg * H0_deg Hadd_deg = H0_deg # parallel to H0 → Schur complement = 0 → A = 0 Wadd_deg, a_deg = GsvdInitialization.init_W(X_deg, W0_deg, H0_deg, Hadd_deg) @test all(isfinite, Wadd_deg) @@ -151,39 +153,39 @@ end U = rand(rng, 10, 3) V = rand(rng, 3, 8) Xdense = U * V - Xfact = MockFactored(U, V) + Xfact = MockFactored(U, V) W0, H0 = rand(rng, 10, 4), rand(rng, 4, 8) - Hadd = rand(rng, 2, 8) - fs = svd(Xdense) - f = (fs.U, fs.S, fs.V) + Hadd = rand(rng, 2, 8) + fs = svd(Xdense) + f = (fs.U, fs.S, fs.V) # `init_W` agrees across both X representations. Wadd_d, a_d = GsvdInitialization.init_W(Xdense, W0, H0, Hadd) - Wadd_f, a_f = GsvdInitialization.init_W(Xfact, W0, H0, Hadd) + Wadd_f, a_f = GsvdInitialization.init_W(Xfact, W0, H0, Hadd) @test Wadd_d ≈ Wadd_f - @test a_d ≈ a_f + @test a_d ≈ a_f # `Wcols_modification` likewise. β0 = rand(rng, 4) W_scaled = repeat(β0', size(W0, 1)) .* W0 @test GsvdInitialization.Wcols_modification(Xdense, W_scaled, H0) ≈ - GsvdInitialization.Wcols_modification(Xfact, W_scaled, H0) + GsvdInitialization.Wcols_modification(Xfact, W_scaled, H0) # End-to-end `gsvdrecover` agrees on the components it returns. Wd, Hd, _ = GsvdInitialization.gsvdrecover(Xdense, copy(W0), copy(H0), 2, f) - Wf, Hf, _ = GsvdInitialization.gsvdrecover(Xfact, copy(W0), copy(H0), 2, f) + Wf, Hf, _ = GsvdInitialization.gsvdrecover(Xfact, copy(W0), copy(H0), 2, f) @test Wd ≈ Wf @test Hd ≈ Hf # The `joint_nnls` strategy also accepts a factored `X` (it needs `X - W*H` # and `eltype(X)` on top of `*`/`sum(abs2, ·)`). Wjd, ajd = GsvdInitialization.init_W_joint_nnls(Xdense, W0, H0, Hadd) - Wjf, ajf = GsvdInitialization.init_W_joint_nnls(Xfact, W0, H0, Hadd) + Wjf, ajf = GsvdInitialization.init_W_joint_nnls(Xfact, W0, H0, Hadd) @test Wjd ≈ Wjf @test ajd ≈ ajf Wd_j, Hd_j, _ = GsvdInitialization.gsvdrecover(GsvdInitialization.joint_nnls, Xdense, copy(W0), copy(H0), 2, f) - Wf_j, Hf_j, _ = GsvdInitialization.gsvdrecover(GsvdInitialization.joint_nnls, Xfact, copy(W0), copy(H0), 2, f) + Wf_j, Hf_j, _ = GsvdInitialization.gsvdrecover(GsvdInitialization.joint_nnls, Xfact, copy(W0), copy(H0), 2, f) @test Wd_j ≈ Wf_j @test Hd_j ≈ Hf_j end @@ -208,14 +210,14 @@ end # as plain `Int` — no `MethodError` from over-tight signatures. Each pair # runs the same pipeline twice, so `rtol = 1e-6` for the same reason as the # cross-call checks in "test top wrapper". - r_int, _ = gsvdnmf(X, 9 => 10; alg = :cd) + r_int, _ = gsvdnmf(X, 9 => 10; alg = :cd) r_int32, _ = gsvdnmf(X, Int32(9) => Int32(10); alg = :cd) - @test isapprox(r_int.W, r_int32.W; rtol = 1e-6) - @test isapprox(r_int.H, r_int32.H; rtol = 1e-6) + @test isapprox(r_int.W, r_int32.W; rtol = 1.0e-6) + @test isapprox(r_int.H, r_int32.H; rtol = 1.0e-6) - r_n2_int, _ = gsvdnmf(X, 10; alg = :cd) + r_n2_int, _ = gsvdnmf(X, 10; alg = :cd) r_n2_int32, _ = gsvdnmf(X, Int32(10); alg = :cd) - @test isapprox(r_n2_int.W, r_n2_int32.W; rtol = 1e-6) + @test isapprox(r_n2_int.W, r_n2_int32.W; rtol = 1.0e-6) # `gsvdrecover` likewise accepts a non-`Int` `kadd`. Wfit, Hfit = rand(rng, 30, 4), rand(rng, 4, 20) @@ -232,15 +234,15 @@ end # BLAS in `nnmf`; `rtol = 1e-6` stays far tighter than the `1e-4` NMF tol. r_default, _ = gsvdnmf(X, 9 => 10; alg = :cd) r_explicit, _ = gsvdnmf(GsvdInitialization.truncating, X, 9 => 10; alg = :cd) - @test isapprox(r_default.W, r_explicit.W; rtol = 1e-6) - @test isapprox(r_default.H, r_explicit.H; rtol = 1e-6) + @test isapprox(r_default.W, r_explicit.W; rtol = 1.0e-6) + @test isapprox(r_default.H, r_explicit.H; rtol = 1.0e-6) # do-block form: anonymous strategy that simply forwards to `truncating` r_doblock, _ = gsvdnmf(X, 9 => 10; alg = :cd) do X0, W0, H0, Hadd GsvdInitialization.truncating(X0, W0, H0, Hadd) end - @test isapprox(r_doblock.W, r_default.W; rtol = 1e-6) - @test isapprox(r_doblock.H, r_default.H; rtol = 1e-6) + @test isapprox(r_doblock.W, r_default.W; rtol = 1.0e-6) + @test isapprox(r_doblock.H, r_default.H; rtol = 1.0e-6) end @testset "integer-n2 convenience methods and :multmse" begin @@ -253,8 +255,8 @@ end r_default, _ = gsvdnmf(X, W, H, 5; alg = :cd) r_strategy, _ = gsvdnmf(GsvdInitialization.truncating, X, W, H, 5; alg = :cd) @test size(r_default.W, 2) == 5 - @test isapprox(r_default.W, r_strategy.W; rtol = 1e-6) - @test isapprox(r_default.H, r_strategy.H; rtol = 1e-6) + @test isapprox(r_default.W, r_strategy.W; rtol = 1.0e-6) + @test isapprox(r_default.H, r_strategy.H; rtol = 1.0e-6) # :multmse floors the augmented factors to `truncmult` so multiplicative # updates (which cannot move entries off zero) can polish them. @@ -316,37 +318,37 @@ end # ~1e-9 nondeterminism of multithreaded BLAS; `rtol = 1e-6` as in the # cross-call checks above. r_plain, _ = gsvdnmf(X, copy(W), copy(H), f; n2 = 5, alg = :cd) - r_view, _ = gsvdnmf(vX, vW, vH, f; n2 = 5, alg = :cd) - @test isapprox(r_plain.W, r_view.W; rtol = 1e-6) - @test isapprox(r_plain.H, r_view.H; rtol = 1e-6) + r_view, _ = gsvdnmf(vX, vW, vH, f; n2 = 5, alg = :cd) + @test isapprox(r_plain.W, r_view.W; rtol = 1.0e-6) + @test isapprox(r_plain.H, r_view.H; rtol = 1.0e-6) end @testset "joint optimize W and alpha" begin rng = StableRNG(7) W = W_GT H = H_GT - X = W*H - result_joint, _ = gsvdnmf(GsvdInitialization.joint_nnls, X, 9=>10; alg = :cd, maxiter = 10^5, tol_final=1e-4, tol_intermediate = 1e-4); + X = W * H + result_joint, _ = gsvdnmf(GsvdInitialization.joint_nnls, X, 9 => 10; alg = :cd, maxiter = 10^5, tol_final = 1.0e-4, tol_intermediate = 1.0e-4) W_gsvd, H_gsvd = result_joint.W, result_joint.H @test size(W_gsvd, 2) == 10 - @test sum(abs2, X-W_gsvd*H_gsvd)/sum(abs2, X) < 2e-10 + @test sum(abs2, X - W_gsvd * H_gsvd) / sum(abs2, X) < 2.0e-10 W, H = rand(rng, 10, 3), rand(rng, 3, 8) - X = W*H + X = W * H U, S, V = svd(X) W0, H0 = copy(W), copy(H) Hadd = rand(rng, 2, 8) Wadd, a = GsvdInitialization.init_W_joint_nnls(X, W0, H0, Hadd) @test a ≈ ones(size(W0, 2)) - @test norm(Wadd) <= 1e-8 + @test norm(Wadd) <= 1.0e-8 G = GsvdInitialization.gram_sp_C(W0, H0, Hadd)[1] b = GsvdInitialization.gram_b(X, W0, H0, Hadd) Wadd = rand(rng, 10, 2) α = rand(rng, 3) θ = vcat(vec(Wadd), α) - E = θ'*G*θ-2*b'*θ+sum(abs2, X) - @test abs(E-sum(abs2, X-[repeat(α', size(W0, 1)).*W0 Wadd]*[H0;Hadd])) <= 1e-12 + E = θ' * G * θ - 2 * b' * θ + sum(abs2, X) + @test abs(E - sum(abs2, X - [repeat(α', size(W0, 1)) .* W0 Wadd] * [H0;Hadd])) <= 1.0e-12 end