Merge branch 'jg/krylovkit'

Project.toml

@@ -5,6 +5,7 @@ version = "1.4.1"
 
 [deps]
 Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
+KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -14,9 +15,10 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
 julia = "1"
+Arpack = "0.3.2"
+KrylovKit = "0.5.2"
 MathOptInterface = "0.9.14"
 TimerOutputs = "0.5.0"
-Arpack = "0.3.2"
 
 [extras]
 

src/ProxSDP.jl

@@ -1,13 +1,15 @@
 module ProxSDP
 
+    using Arpack
+    using KrylovKit
     using MathOptInterface
     using TimerOutputs
-    using Arpack
+
     using Printf
     using SparseArrays
     using LinearAlgebra
-
     import Random
+
     import LinearAlgebra: BlasInt
 
     include("structs.jl")

src/eigsolver.jl

@@ -729,22 +729,47 @@ end
 function eig!(arc::ARPACKAlloc, A::Symmetric{T,Matrix{T}}, nev::Integer, opt::Options)::Nothing where T
 
     up_ncv = false
-    for i in 1:1
-        # Initialize parameters and do memory allocation
-        @timeit "update_arc" _update_arc!(arc, A, nev, opt, up_ncv)::Nothing
+    if opt.eigsolver == 1
+        for i in 1:1
+            # Initialize parameters and do memory allocation
+            @timeit "update_arc" _update_arc!(arc, A, nev, opt, up_ncv)::Nothing
+
+            # Top level reverse communication interface to solve real double precision symmetric problems.
+            @timeit "saupd" _saupd!(arc)::Nothing
+
+            if arc.nev > arc.iparam[5]
+                up_ncv = true
+            else
+                break
+            end
+        end
 
-        # Top level reverse communication interface to solve real double precision symmetric problems.
-        @timeit "saupd" _saupd!(arc)::Nothing
+        # Post processing routine (eigenvalues and eigenvectors purification)
+        @timeit "seupd" _seupd!(arc)::Nothing
+    else
+        @timeit "update_arc" _update_arc!(arc, A, nev, opt, up_ncv)::Nothing
+        @timeit "krylovkit" begin
+            vals, vecs, info = eigsolve(
+                A, arc.resid, arc.nev, :LR,Lanczos(krylovdim = arc.ncv, tol = 1e-12, maxiter=100)
+            )
 
-        if arc.nev > arc.iparam[5]
-            up_ncv = true
-        else
-            break
+            if info.converged == 0
+                @show arc.nev, info
+            end
+            fill!(arc.d, 0.0)
+            fill!(arc.v, 0.0)
+            for i in 1:min(arc.nev, info.converged)
+                arc.d[i] = vals[i]
+                for j in 1:arc.n
+                    arc.v[j,i] = vecs[i][j]
+                end
+            end
+            arc.converged = true
+            arc.arpackerror = false
+            # if info.converged == 1
+            #     arc.converged = false
+            # end
         end
     end
-
-    # Post processing routine (eigenvalues and eigenvectors purification)
-    @timeit "seupd" _seupd!(arc)::Nothing
-
     return nothing
 end

src/prox_operators.jl

@@ -27,15 +27,19 @@ function psd_projection!(v::Vector{Float64}, a::AuxiliaryData, cones::ConicSets,
             a.m[idx][1] = max(0., a.m[idx][1])
             p.min_eig[idx] = a.m[idx][1]
 
-        elseif !opt.full_eig_decomp && p.target_rank[idx] <= opt.max_target_rank_krylov_eigs && sdp.sq_side > opt.min_size_krylov_eigs
+        elseif !opt.full_eig_decomp &&
+            p.target_rank[idx] <= opt.max_target_rank_krylov_eigs &&
+            sdp.sq_side > opt.min_size_krylov_eigs &&
+            mod(p.iter, opt.full_eig_freq) > opt.full_eig_len # for full from time to time
             @timeit "eigs" begin 
                 eig!(arc_list[idx], a.m[idx], p.target_rank[idx], opt)
+
                 if hasconverged(arc_list[idx])
                     fill!(a.m[idx].data, 0.)
                     for i in 1:p.target_rank[idx]
                         if unsafe_getvalues(arc_list[idx])[i] > 0. 
                             p.current_rank[idx] += 1
-                            vec = unsafe_getvectors(arc_list[idx])[:, i]
+                            vec = view(unsafe_getvectors(arc_list[idx]), :, i)
                             LinearAlgebra.BLAS.gemm!('N', 'T', unsafe_getvalues(arc_list[idx])[i], vec, vec, 1., a.m[idx].data)
                         end
                     end
@@ -74,7 +78,8 @@ function full_eig!(a::AuxiliaryData, idx::Int, opt::Options, p::Params)::Nothing
     fill!(a.m[idx].data, 0.)
     for i in 1:length(fact.values)
         if fact.values[i] > 0.
-            LinearAlgebra.BLAS.gemm!('N', 'T', fact.values[i], fact.vectors[:, i], fact.vectors[:, i], 1., a.m[idx].data)
+            v = view(fact.vectors, :, i)
+            LinearAlgebra.BLAS.gemm!('N', 'T', fact.values[i], v, v, 1., a.m[idx].data)
             if fact.values[i] > opt.tol_psd
                 p.current_rank[idx] += 1
             end

src/structs.jl

@@ -51,9 +51,9 @@ Base.@kwdef mutable struct Options
 
     # Linesearch parameters
     max_linsearch_steps::Int = 5000
-    delta::Float64 = .999
+    delta::Float64 = .9999
     initial_theta::Float64 = 1.
-    linsearch_decay::Float64 = .9
+    linsearch_decay::Float64 = .75
 
     # Spectral decomposition parameters
     full_eig_decomp::Bool = false
@@ -63,6 +63,13 @@ Base.@kwdef mutable struct Options
     rank_increment::Int = 1 # 0=multiply, 1 = add
     rank_increment_factor::Int = 1 # 0 multiply, 1 = add
 
+    # eigsolver selection
+    #=
+        1: Arpack [dsaupd] (tipically non-deterministic)
+        2: KrylovKit [eigsolve/Lanczos] (DEFAULT)
+    =#
+    eigsolver::Int = 2
+
     # Arpack
     arpack_tol::Float64 = 1e-10
     #=
@@ -73,7 +80,7 @@ Base.@kwdef mutable struct Options
     =#
     arpack_resid_init::Int = 3
     arpack_resid_seed::Int = 1234
-    arpack_reset_resid::Bool = true # true for determinism
+    arpack_reset_resid::Bool = false # true for determinism
     arpack_max_iter::Int = 10_000
     arpack_min_lanczos::Int = 25
     # larger is more stable to converge and more deterministic
@@ -83,6 +90,9 @@ Base.@kwdef mutable struct Options
     reduce_rank::Bool = false
     rank_slack::Int = 3
 
+    full_eig_freq::Int = 10000000
+    full_eig_len::Int = 0
+
     # equilibration parameters
     equilibration::Bool = true
     equilibration_iters::Int = 1000