Some normalized_cut() fixes (#84)

* eigs(): cleanup, sort complex eigvals by real part

* normalized_cut(): check size of weightmatrix early

* normalized_cut(): take care of isolated vertices

because they result in degenerated matrix

* normalized_cut(): check if 2nd eigvec is missing

* normalized_cut(): use real part of eigenvec

---------

Co-authored-by: Guillaume Dalle <[email protected]>

Author: alyst

src/graphcut/normalized_cut.jl

@@ -117,26 +117,39 @@ end
 
 function _recursive_normalized_cut(W, thres, num_cuts)
     m, n = size(W)
-    D = Diagonal(vec(sum(W; dims=2)))
-
-    m == 1 && return [1]
+    (m <= 1) && return ones(Int, m) # trivial
+    D = Diagonal(vec(sum(W, dims=2)))
+
+    # check that the diagonal is not degenerated as otherwise invDroot errors
+    dnz = abs.(diag(D)) .>= 1E-16
+    if !all(dnz)
+        # vertices with incident edges summing to almost zero
+        # are not connected to the rest of the subnetwork,
+        # put them to separate modules and cut the remaining submatrix
+        nzlabels = _recursive_normalized_cut(W[dnz, dnz], thres, num_cuts)
+        nzix = 0
+        zix = maximum(nzlabels)
+        return Int[nz ? nzlabels[nzix += 1] : (zix += 1) for nz in dnz]
+    end
 
-    # get eigenvector corresponding to second smallest eigenvalue
+    # get eigenvector corresponding to the second smallest generalized eigenvalue:
     # v = eigs(D-W, D, nev=2, which=SR())[2][:,2]
     # At least some versions of ARPACK have a bug, this is a workaround
     invDroot = sqrt.(inv(D)) # equal to Cholesky factorization for diagonal D
     if n > 12
-        λ, Q = eigs(invDroot' * (D - W) * invDroot; nev=12, which=SR())
-        ret = real(Q[:, 2])
+        _, Q = eigs(invDroot' * (D - W) * invDroot, nev=12, which=SR())
+        (size(Q, 2) <= 1) && return collect(1:m) # no 2nd eigenvector
+        ret = convert(Vector, real(view(Q, :, 2)))
     else
         ret = eigen(Matrix(invDroot' * (D - W) * invDroot)).vectors[:, 2]
     end
-    v = invDroot * ret
+    v = real(invDroot * ret)
 
     # perform n-cuts with different partitions of v and find best one
     min_cost = Inf
     best_thres = -1
-    for t in range(minimum(v); stop=maximum(v), length=num_cuts)
+    vmin, vmax = extrema(v)
+    for t in range(vmin, stop=vmax, length=num_cuts)
         cut = v .> t
         cost = _normalized_cut_cost(cut, W, D)
         if cost < min_cost

src/linalg/LinAlg.jl

@@ -51,13 +51,8 @@ function eigs(A; kwargs...)
     schr = partialschur(A; kwargs...)
     vals, vectors = partialeigen(schr[1])
     reved = (kwargs[:which] == LR() || kwargs[:which] == LM())
-    k::Int = get(kwargs, :nev, length(vals))
-    k = min(k, length(vals))
-    perm = collect(1:k)
-    if vals[1] isa (Real)
-        perm = sortperm(vals; rev=reved)
-        perm = perm[1:k]
-    end
+    k = min(get(kwargs, :nev, length(vals))::Int, length(vals))
+    perm = sortperm(vals, by=real, rev=reved)[1:k]
     λ = vals[perm]
     Q = vectors[:, perm]
     return λ, Q