deal with Inf case, change jump (#2)

* deal with Inf case, change jump

* Update src/GsvdInitialization.jl

Co-authored-by: Tim Holy <tim.holy@gmail.com>

---------

Co-authored-by: youdongguo <1010705897@qq.com>
Co-authored-by: Tim Holy <tim.holy@gmail.com>

Author: 3 people

Project.toml

@@ -4,14 +4,11 @@ authors = ["youdongguo <1010705897@qq.com> and contributors"]
 version = "0.1.0"
 
 [deps]
-Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
-JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NMF = "6ef6ca0d-6ad7-5ff6-b225-e928bfa0a386"
+NonNegLeastSquares = "b7351bd1-99d9-5c5d-8786-f205a815c4d7"
 
 [compat]
-Ipopt = "1"
-JuMP = "1"
 NMF = "1"
 julia = "1"
 

src/GsvdInitialization.jl

@@ -1,74 +1,70 @@
 module GsvdInitialization
 
 using LinearAlgebra, NMF
-using JuMP, Ipopt
+using NonNegLeastSquares
 
-export overnmfinit,
-       gsvdinit,
-       init_H,
-       init_W,
-       Wcols_modification
+export gsvdnmf,
+       gsvdrecover
 
-function overnmfinit(X::AbstractArray, W0::AbstractArray, H0::AbstractArray, kadd::Int; initdata = nothing, n::Int = size(W0, 2))
+function gsvdnmf(X, ncomponents::Pair{Int,Int}; tol_final=1e-4, tol_intermediate=sqrt(tol_final), W0=nothing, H0=nothing, kwargs...)
+    f = svd(X)
+    if W0 === nothing && H0 === nothing
+        W0, H0 = NMF.nndsvd(X, ncomponents[2], initdata=f)
+    end
+    result_initial = nnmf(X, ncomponents[2]; kwargs..., init=:custom, tol=tol_intermediate, W0=copy(W0), H0=copy(H0))
+    W_initial, H_initial = result_initial.W, result_initial.H
+    kadd = ncomponents[2] - ncomponents[1]
+    kadd >= 0 || throw(ArgumentError("The number of components to add must be non-negative."))
+    kadd <= ncomponents[2] || throw(ArgumentError("The number of components to add must be less than the total number of components."))
+    W_recover, H_recover = gsvdrecover(X, copy(W_initial), copy(H_initial), kadd, initdata=f)
+    result_recover = nnmf(X, ncomponents[1]; kwargs..., init=:custom, tol=tol_final, W0=copy(W_recover), H0=copy(H_recover))
+    return result_recover
+end
+    
+function gsvdrecover(X::AbstractArray, W0::AbstractArray, H0::AbstractArray, kadd::Int; initdata = nothing)
     if kadd == 0
         return W0, H0
     else
         m = size(W0, 1) 
-        Wadd, Hadd, a = gsvdinit(X, W0, H0, kadd; initdata = initdata, n = n)
+        Wadd, Hadd, a, Λ = components_recover(X, W0, H0, kadd; initdata = initdata)
         Wadd_nn, Hadd_nn = NMF.nndsvd(X, kadd, initdata = (U = Wadd, S = ones(kadd), V = Hadd'))
         W0_1, H0_1 = [repeat(a', m, 1).*W0 Wadd_nn], [H0; Hadd_nn]
         cs = Wcols_modification(X, W0_1, H0_1)
         W0_2, H0_2 = repeat(cs', m, 1).*W0_1, H0_1
-        return abs.(W0_2), abs.(H0_2)
+        return abs.(W0_2), abs.(H0_2), Λ
     end
 end
 
-function gsvdinit(X::AbstractArray, W0::AbstractArray, H0::AbstractArray, kadd::Int; initdata = nothing, n::Int = size(W0, 2))
-        # n = size(W0, 2)+kadd
-        # @show n
+function components_recover(X::AbstractArray, W0::AbstractArray, H0::AbstractArray, kadd::Int; initdata = nothing)
+        n::Int = size(W0, 2)
         kadd <= n || throw(ArgumentError("# of extra columns must less than 1st NMF components"))
         U, S, V = initdata === nothing ? svd(X) : (initdata.U, initdata.S, initdata.V)
         U0, S0, V0 = U[:,1:n], S[1:n], V[:,1:n]
-        Hadd = init_H(U0, S0, V0, W0, H0, kadd)
+        Hadd, Λ = init_H(U0, S0, V0, W0, H0, kadd)
         Wadd, a = init_W(X, W0, H0, Hadd)
-        return Wadd, Hadd, a
+        return Wadd, Hadd, a, Λ
 end
 
 function init_H(U0::AbstractArray, S0::AbstractArray, V0::AbstractArray, W0::AbstractArray, H0::AbstractArray, kadd::Int)
     _, _, Q, D1, D2, R = svd(Matrix(Diagonal(S0)), (U0'*W0)*(H0*V0));
-    # _, _, Q, D1, D2, R = svd(Matrix(Diagonal(S0)), W0*(H0*V0));
     inv_RQt = inv(R*Q')
-    HHH = inv_RQt
-    F = (diag(D1)./diag(D2)).^2
-    # @show F
-    if kadd < size(U0, 2)
-        k0 = kadd
-        H_index = Int[]
-        while k0 >= 1
-            j = findmax(F)[2]
-            F[j] = -1
-            push!(H_index, j) 
-            k0 -= 1   
-        end
-        Hadd = HHH[:,H_index]
-    else
-        Hadd = HHH
-    end
+    r0 = size(U0, 2)
+    k = findfirst(x->x!=0, D2[1,:])
+    k = (k === nothing) ? r0 : k-1
+    kadd >= k || @warn "kadd is less than rank deficiency of W0*H0."
+    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]
+    Hadd = inv_RQt[:, H_index]
     Hadd_1 = V0*Hadd
-    return Hadd_1'
+    return Hadd_1', Λ[H_index]
 end
 
 function init_W(X::AbstractArray{T}, W0::AbstractArray{T}, H0::AbstractArray{T}, Hadd::AbstractArray{T}; α = nothing) where T
-    R = size(W0, 2)
     A, b, _, invHH, H0Hadd, XHaddt = obj_para(X, W0, H0, Hadd)
-    if α === nothing 
-        model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
-        @variable(model, a[1:R] >= 1e-12, start = 1)
-        @objective(model, Min, a'*A*a+2*b'*a)
-        optimize!(model)
-        α = JuMP.value.(a)
-    end
-    Wadd = XHaddt*invHH-W0*Diagonal(α)*H0Hadd*invHH
+    (isposdef(A) || sum(abs2, A) <= 1e-12) || @warn "A is not positive definite."
+    α = α === nothing ? nonneg_lsq(A, -b; alg=:fnnls, gram=true) : α
+    Wadd = XHaddt*invHH-W0*Diagonal(α[:])*H0Hadd*invHH
     return Wadd, abs.(α)
 end
 
@@ -84,7 +80,7 @@ function obj_para(X::AbstractArray{T}, W0::AbstractArray{T}, H0::AbstractArray{T
     W0XHaddt = W0'*XHaddt
     b = diag(H0Hadd*invHH*W0XHaddt'-W0tXH0t)
     C = sum(abs2, X)-sum(invHH.*(XHaddt'*XHaddt))
-    return A, b, C, invHH, H0Hadd, XHaddt
+    return Symmetric(A), b, C, invHH, H0Hadd, XHaddt
 end
 
 function Wcols_modification(X::AbstractArray{T}, W::AbstractArray{T}, H::AbstractArray{T}) where T
@@ -95,12 +91,8 @@ function Wcols_modification(X::AbstractArray{T}, W::AbstractArray{T}, H::Abstrac
     WtXHt = W'*X*H'
     a = diag(WtXHt)
     B = WW.*HH
-    model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
-    @variable(model, b[1:n] >= 1e-12, start = 1)
-    @objective(model, Min, b'*B*b-2*a'*b)
-    optimize!(model)
-    β = JuMP.value.(b)
-    return β
+    β = nonneg_lsq(B, a; alg=:fnnls, gram=true) 
+    return β[:]
 end
 
 end

test/runtests.jl

@@ -24,12 +24,12 @@ using LinearAlgebra, NMF
     Hadd = rand(2, 8)
     A, b, C, HH, γ = GsvdInitialization.obj_para(X, W0, H0, Hadd)
     a = rand(4)
-    Wadd, a = init_W(X, W0, H0, Hadd, α = a)
+    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
 
     β0 = rand(3)
-    β = Wcols_modification(X, repeat(β0', size(W, 1)).*W, H)
+    β = GsvdInitialization.Wcols_modification(X, repeat(β0', size(W, 1)).*W, H)
     @test β.*β0 ≈ ones(3)
 
 end