engine-proto/gram-test/Engine.jl

@@ -1,14 +1,10 @@
 module Engine
 
 using LinearAlgebra
-using GenericLinearAlgebra
 using SparseArrays
 using Random
-using Optim
 
-export
+export rand_on_shell, Q, DescentHistory, realize_gram
-  rand_on_shell, Q, DescentHistory,
-  realize_gram_gradient, realize_gram_newton, realize_gram_optim, realize_gram
 
 # === guessing ===
 
@@ -80,11 +76,8 @@ end
 struct DescentHistory{T}
   scaled_loss::Array{T}
   neg_grad::Array{Matrix{T}}
-  base_step::Array{Matrix{T}}
-  hess::Array{Hermitian{T, Matrix{T}}}
   slope::Array{T}
   stepsize::Array{T}
-  positive::Array{Bool}
   backoff_steps::Array{Int64}
   last_line_L::Array{Matrix{T}}
   last_line_loss::Array{T}
@@ -92,16 +85,13 @@ struct DescentHistory{T}
   function DescentHistory{T}(
     scaled_loss = Array{T}(undef, 0),
     neg_grad = Array{Matrix{T}}(undef, 0),
-    hess = Array{Hermitian{T, Matrix{T}}}(undef, 0),
-    base_step = Array{Matrix{T}}(undef, 0),
     slope = Array{T}(undef, 0),
     stepsize = Array{T}(undef, 0),
-    positive = Bool[],
     backoff_steps = Int64[],
     last_line_L = Array{Matrix{T}}(undef, 0),
     last_line_loss = Array{T}(undef, 0)
   ) where T
-    new(scaled_loss, neg_grad, hess, base_step, slope, stepsize, positive, backoff_steps, last_line_L, last_line_loss)
+    new(scaled_loss, neg_grad, slope, stepsize, backoff_steps, last_line_L, last_line_loss)
   end
 end
 
@@ -111,7 +101,7 @@ function realize_gram_gradient(
   gram::SparseMatrixCSC{T, <:Any},
   guess::Matrix{T};
   scaled_tol = 1e-30,
-  min_efficiency = 0.5,
+  target_improvement = 0.5,
   init_stepsize = 1.0,
   backoff = 0.9,
   max_descent_steps = 600,
@@ -162,7 +152,7 @@ function realize_gram_gradient(
       improvement = loss_last - loss
       push!(history.last_line_L, L)
       push!(history.last_line_loss, loss / scale_adjustment)
-      if improvement >= min_efficiency * stepsize * slope
+      if improvement >= target_improvement * stepsize * slope
         history.backoff_steps[end] = backoff_steps
         break
       end
@@ -211,7 +201,7 @@ function realize_gram_newton(
   scale_adjustment = sqrt(T(length(constrained)))
   tol = scale_adjustment * scaled_tol
   
-  # use Newton's method
+  # use newton's method
   L = copy(guess)
   for step in 0:max_steps
     # evaluate the loss function
@@ -239,10 +229,8 @@ function realize_gram_newton(
       deriv_grad = 4*Q*(-basis_mat*Δ_proj + L*neg_dΔ_proj)
       hess[:, (k-1)*element_dim + j] = reshape(deriv_grad, total_dim)
     end
-    hess = Hermitian(hess)
-    push!(history.hess, hess)
     
-    # compute the Newton step
+    # compute the newton step
     step = hess \ reshape(neg_grad, total_dim)
     L += rate * reshape(step, dims)
   end
@@ -251,168 +239,4 @@ function realize_gram_newton(
   L, history
 end
 
-LinearAlgebra.eigen!(A::Symmetric{BigFloat, Matrix{BigFloat}}; sortby::Nothing) =
-  eigen!(Hermitian(A))
-
-function convertnz(type, mat)
-  J, K, values = findnz(mat)
-  sparse(J, K, type.(values))
-end
-
-function realize_gram_optim(
-  gram::SparseMatrixCSC{T, <:Any},
-  guess::Matrix{T}
-) where T <: Number
-  # find the dimension of the search space
-  dims = size(guess)
-  element_dim, construction_dim = dims
-  total_dim = element_dim * construction_dim
-  
-  # list the constrained entries of the gram matrix
-  J, K, _ = findnz(gram)
-  constrained = zip(J, K)
-  
-  # scale the loss function
-  scale_adjustment = length(constrained)
-  
-  function loss(L_vec)
-    L = reshape(L_vec, dims)
-    Δ_proj = proj_diff(gram, L'*Q*L)
-    dot(Δ_proj, Δ_proj) / scale_adjustment
-  end
-  
-  function loss_grad!(storage, L_vec)
-    L = reshape(L_vec, dims)
-    Δ_proj = proj_diff(gram, L'*Q*L)
-    storage .= reshape(-4*Q*L*Δ_proj, total_dim) / scale_adjustment
-  end
-  
-  function loss_hess!(storage, L_vec)
-    L = reshape(L_vec, dims)
-    Δ_proj = proj_diff(gram, L'*Q*L)
-    indices = [(j, k) for k in 1:construction_dim for j in 1:element_dim]
-    for (j, k) in indices
-      basis_mat = basis_matrix(T, j, k, dims)
-      neg_dΔ = basis_mat'*Q*L + L'*Q*basis_mat
-      neg_dΔ_proj = proj_to_entries(neg_dΔ, constrained)
-      deriv_grad = 4*Q*(-basis_mat*Δ_proj + L*neg_dΔ_proj) / scale_adjustment
-      storage[:, (k-1)*element_dim + j] = reshape(deriv_grad, total_dim)
-    end
-  end
-  
-  optimize(
-    loss, loss_grad!, loss_hess!,
-    reshape(guess, total_dim),
-    Newton()
-  )
-end
-
-# seek a matrix `L` for which `L'QL` matches the sparse matrix `gram` at every
-# explicit entry of `gram`. use gradient descent starting from `guess`
-function realize_gram(
-  gram::SparseMatrixCSC{T, <:Any},
-  guess::Matrix{T};
-  scaled_tol = 1e-30,
-  min_efficiency = 0.5,
-  init_rate = 1.0,
-  backoff = 0.9,
-  reg_scale = 1.1,
-  max_descent_steps = 200,
-  max_backoff_steps = 110
-) where T <: Number
-  # start history
-  history = DescentHistory{T}()
-  
-  # find the dimension of the search space
-  dims = size(guess)
-  element_dim, construction_dim = dims
-  total_dim = element_dim * construction_dim
-  
-  # list the constrained entries of the gram matrix
-  J, K, _ = findnz(gram)
-  constrained = zip(J, K)
-  
-  # scale the tolerance
-  scale_adjustment = sqrt(T(length(constrained)))
-  tol = scale_adjustment * scaled_tol
-  
-  # initialize variables
-  grad_rate = init_rate
-  L = copy(guess)
-  
-  # use Newton's method with backtracking and gradient descent backup
-  Δ_proj = proj_diff(gram, L'*Q*L)
-  loss = dot(Δ_proj, Δ_proj)
-  for step in 1:max_descent_steps
-    # stop if the loss is tolerably low
-    if loss < tol
-      break
-    end
-    
-    # find the negative gradient of loss function
-    neg_grad = 4*Q*L*Δ_proj
-    
-    # find the negative Hessian of the loss function
-    hess = Matrix{T}(undef, total_dim, total_dim)
-    indices = [(j, k) for k in 1:construction_dim for j in 1:element_dim]
-    for (j, k) in indices
-      basis_mat = basis_matrix(T, j, k, dims)
-      neg_dΔ = basis_mat'*Q*L + L'*Q*basis_mat
-      neg_dΔ_proj = proj_to_entries(neg_dΔ, constrained)
-      deriv_grad = 4*Q*(-basis_mat*Δ_proj + L*neg_dΔ_proj)
-      hess[:, (k-1)*element_dim + j] = reshape(deriv_grad, total_dim)
-    end
-    hess = Hermitian(hess)
-    push!(history.hess, hess)
-    
-    # regularize the Hessian
-    min_eigval = minimum(eigvals(hess))
-    push!(history.positive, min_eigval > 0)
-    if min_eigval <= 0
-      hess -= reg_scale * min_eigval * I
-    end
-    base_step = reshape(hess \ reshape(neg_grad, total_dim), dims)
-    push!(history.base_step, base_step)
-    
-    # store the current position, loss, and slope
-    L_last = L
-    loss_last = loss
-    push!(history.scaled_loss, loss / scale_adjustment)
-    push!(history.neg_grad, neg_grad)
-    push!(history.slope, norm(neg_grad))
-    
-    # find a good step size using backtracking line search
-    push!(history.stepsize, 0)
-    push!(history.backoff_steps, max_backoff_steps)
-    empty!(history.last_line_L)
-    empty!(history.last_line_loss)
-    rate = one(T)
-    step_success = false
-    for backoff_steps in 0:max_backoff_steps
-      history.stepsize[end] = rate
-      L = L_last + rate * base_step
-      Δ_proj = proj_diff(gram, L'*Q*L)
-      loss = dot(Δ_proj, Δ_proj)
-      improvement = loss_last - loss
-      push!(history.last_line_L, L)
-      push!(history.last_line_loss, loss / scale_adjustment)
-      if improvement >= min_efficiency * rate * dot(neg_grad, base_step)
-        history.backoff_steps[end] = backoff_steps
-        step_success = true
-        break
-      end
-      rate *= backoff
-    end
-    
-    # if we've hit a wall, quit
-    if !step_success
-      return L_last, false, history
-    end
-  end
-  
-  # return the factorization and its history
-  push!(history.scaled_loss, loss / scale_adjustment)
-  L, true, history
-end
-
 end

engine-proto/gram-test/circles-in-triangle.jl

@@ -81,16 +81,12 @@ guess = hcat(
 =#
 
 # complete the gram matrix using gradient descent followed by Newton's method
-#=
 L, history = Engine.realize_gram_gradient(gram, guess, scaled_tol = 0.01)
 L_pol, history_pol = Engine.realize_gram_newton(gram, L, rate = 0.3, scaled_tol = 1e-9)
 L_pol2, history_pol2 = Engine.realize_gram_newton(gram, L_pol)
-=#
-L, success, history = Engine.realize_gram(gram, guess)
 completed_gram = L'*Engine.Q*L
 println("Completed Gram matrix:\n")
 display(completed_gram)
-#=
 println(
   "\nSteps: ",
   size(history.scaled_loss, 1),
@@ -98,11 +94,3 @@ println(
   " + ", size(history_pol2.scaled_loss, 1)
 )
 println("Loss: ", history_pol2.scaled_loss[end], "\n")
-=#
-if success
-  println("\nTarget accuracy achieved!")
-else
-  println("\nFailed to reach target accuracy")
-end
-println("Steps: ", size(history.scaled_loss, 1))
-println("Loss: ", history.scaled_loss[end], "\n")

engine-proto/gram-test/overlapping-pyramids.jl

@@ -47,30 +47,17 @@ guess = hcat(
   Engine.rand_on_shell(fill(BigFloat(-1), 2))
 )
 
-# complete the gram matrix
+# complete the gram matrix using gradient descent followed by Newton's method
-#=
 L, history = Engine.realize_gram_gradient(gram, guess, scaled_tol = 0.01)
 L_pol, history_pol = Engine.realize_gram_newton(gram, L)
-=#
-L, success, history = Engine.realize_gram(gram, guess)
 completed_gram = L'*Engine.Q*L
 println("Completed Gram matrix:\n")
 display(completed_gram)
-#=
 println("\nSteps: ", size(history.scaled_loss, 1), " + ", size(history_pol.scaled_loss, 1))
 println("Loss: ", history_pol.scaled_loss[end], "\n")
-=#
-if success
-  println("\nTarget accuracy achieved!")
-else
-  println("\nFailed to reach target accuracy")
-end
-println("Steps: ", size(history.scaled_loss, 1))
-println("Loss: ", history.scaled_loss[end], "\n")
 
 # === algebraic check ===
 
-#=
 R, gens = polynomial_ring(AbstractAlgebra.Rationals{BigInt}(), ["x", "t₁", "t₂", "t₃"])
 x = gens[1]
 t = gens[2:4]
@@ -98,4 +85,3 @@ x_constraint = 25//16 * to_univariate(S, evaluate(rank_constraints[1], [2], [ind
 t₂_constraint = 25//16 * to_univariate(S, evaluate(rank_constraints[3], [2], [indep_val]))
 x_vals = PolynomialRoots.roots(x_constraint.coeffs)
 t₂_vals = PolynomialRoots.roots(t₂_constraint.coeffs)
-=#

engine-proto/gram-test/sphere-in-tetrahedron.jl

@@ -33,18 +33,10 @@ guess = sqrt(1/BigFloat(3)) * BigFloat[
   1  1  1  1  1
 ] + 0.2*Engine.rand_on_shell(fill(BigFloat(-1), 5))
 
-# complete the gram matrix
+# complete the gram matrix using Newton's method
-#=
 L, history = Engine.realize_gram_newton(gram, guess)
-=#
-L, success, history = Engine.realize_gram(gram, guess)
 completed_gram = L'*Engine.Q*L
 println("Completed Gram matrix:\n")
 display(completed_gram)
-if success
+println("\nSteps: ", size(history.scaled_loss, 1))
-  println("\nTarget accuracy achieved!")
-else
-  println("\nFailed to reach target accuracy")
-end
-println("Steps: ", size(history.scaled_loss, 1))
 println("Loss: ", history.scaled_loss[end], "\n")