Skip to content

Apply runic formatting #32

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
timholy merged 1 commit into from
Jun 12, 2026
Merged

Apply runic formatting #32

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
102 changes: 52 additions & 50 deletions src/GsvdInitialization.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"))
Expand Down Expand Up @@ -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
Expand All @@ -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...) =
Expand Down Expand Up @@ -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...)

Expand Down Expand Up @@ -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, Λ
Expand Down Expand Up @@ -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
Expand All @@ -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))
Expand All @@ -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)
Expand All @@ -412,4 +414,4 @@ function truncatepos(Y, X, W, H)
end


end
end
Loading