diff --git a/engine-proto/gram-test/Engine.jl b/engine-proto/gram-test/Engine.jl index e8bda9c..e6f0d97 100644 --- a/engine-proto/gram-test/Engine.jl +++ b/engine-proto/gram-test/Engine.jl @@ -51,11 +51,21 @@ end # the Lorentz form Q = diagm([1, 1, 1, 1, -1]) +# project a matrix onto the subspace of matrices whose entries vanish at the +# given indices +function proj_to_entries(mat, indices) + result = zeros(size(mat)) + for (j, k) in indices + result[j, k] = mat[j, k] + end + result +end + # the difference between the matrices `target` and `attempt`, projected onto the # subspace of matrices whose entries vanish at each empty index of `target` function proj_diff(target::SparseMatrixCSC{T, <:Any}, attempt::Matrix{T}) where T J, K, values = findnz(target) - result = zeros(size(target)...) + result = zeros(size(target)) for (j, k, val) in zip(J, K, values) result[j, k] = val - attempt[j, k] end @@ -65,23 +75,29 @@ end # a type for keeping track of gradient descent history struct DescentHistory{T} scaled_loss::Array{T} + neg_grad::Array{Matrix{T}} slope::Array{T} stepsize::Array{T} backoff_steps::Array{Int64} + last_line_L::Array{Matrix{T}} + last_line_loss::Array{T} function DescentHistory{T}( scaled_loss = Array{T}(undef, 0), + neg_grad = Array{Matrix{T}}(undef, 0), slope = Array{T}(undef, 0), stepsize = Array{T}(undef, 0), - backoff_steps = Int64[] + backoff_steps = Int64[], + last_line_L = Array{Matrix{T}}(undef, 0), + last_line_loss = Array{T}(undef, 0) ) where T - new(scaled_loss, slope, stepsize, backoff_steps) + new(scaled_loss, neg_grad, slope, stepsize, backoff_steps, last_line_L, last_line_loss) end 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( +function realize_gram_gradient( gram::SparseMatrixCSC{T, <:Any}, guess::Matrix{T}; scaled_tol = 1e-30, @@ -104,8 +120,8 @@ function realize_gram( # do gradient descent Δ_proj = proj_diff(gram, L'*Q*L) - loss = norm(Δ_proj) - for step in 1:max_descent_steps + loss = dot(Δ_proj, Δ_proj) + for _ in 1:max_descent_steps # stop if the loss is tolerably low if loss < tol break @@ -114,28 +130,39 @@ function realize_gram( # find negative gradient of loss function neg_grad = 4*Q*L*Δ_proj slope = norm(neg_grad) + dir = neg_grad / slope # store 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, slope) # 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) for backoff_steps in 0:max_backoff_steps history.stepsize[end] = stepsize - L = L_last + stepsize * neg_grad + L = L_last + stepsize * dir Δ_proj = proj_diff(gram, L'*Q*L) - loss = norm(Δ_proj) + loss = dot(Δ_proj, Δ_proj) improvement = loss_last - loss + push!(history.last_line_L, L) + push!(history.last_line_loss, loss / scale_adjustment) if improvement >= target_improvement * stepsize * slope history.backoff_steps[end] = backoff_steps break end stepsize *= backoff end + + # [DEBUG] if we've hit a wall, quit + if history.backoff_steps[end] == max_backoff_steps + break + end end # return the factorization and its history @@ -143,4 +170,73 @@ function realize_gram( L, history end +function basis_matrix(::Type{T}, j, k, dims) where T + result = zeros(T, dims) + result[j, k] = one(T) + result +end + +# seek a matrix `L` for which `L'QL` matches the sparse matrix `gram` at every +# explicit entry of `gram`. use Newton's method starting from `guess` +function realize_gram_newton( + gram::SparseMatrixCSC{T, <:Any}, + guess::Matrix{T}; + scaled_tol = 1e-30, + rate = 1, + max_steps = 100 +) 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 + + # use newton's method + L = copy(guess) + for step in 0:max_steps + # evaluate the loss function + Δ_proj = proj_diff(gram, L'*Q*L) + loss = dot(Δ_proj, Δ_proj) + + # store the current loss + push!(history.scaled_loss, loss / scale_adjustment) + + # stop if the loss is tolerably low + if loss < tol || step > max_steps + 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 + + # compute the newton step + step = hess \ reshape(neg_grad, total_dim) + L += rate * reshape(step, dims) + end + + # return the factorization and its history + L, history +end + end \ No newline at end of file diff --git a/engine-proto/gram-test/basin-shapes.jl b/engine-proto/gram-test/basin-shapes.jl new file mode 100644 index 0000000..5c03c01 --- /dev/null +++ b/engine-proto/gram-test/basin-shapes.jl @@ -0,0 +1,99 @@ +include("Engine.jl") + +using LinearAlgebra +using SparseArrays + +function sphere_in_tetrahedron_shape() + # initialize the partial gram matrix for a sphere inscribed in a regular + # tetrahedron + J = Int64[] + K = Int64[] + values = BigFloat[] + for j in 1:5 + for k in 1:5 + push!(J, j) + push!(K, k) + if j == k + push!(values, 1) + elseif (j <= 4 && k <= 4) + push!(values, -1/BigFloat(3)) + else + push!(values, -1) + end + end + end + gram = sparse(J, K, values) + + # plot loss along a slice + loss_lin = [] + loss_sq = [] + mesh = range(0.9, 1.1, 101) + for t in mesh + L = hcat( + Engine.plane(normalize(BigFloat[ 1, 1, 1]), BigFloat(1)), + Engine.plane(normalize(BigFloat[ 1, -1, -1]), BigFloat(1)), + Engine.plane(normalize(BigFloat[-1, 1, -1]), BigFloat(1)), + Engine.plane(normalize(BigFloat[-1, -1, 1]), BigFloat(1)), + Engine.sphere(BigFloat[0, 0, 0], BigFloat(t)) + ) + Δ_proj = Engine.proj_diff(gram, L'*Engine.Q*L) + push!(loss_lin, norm(Δ_proj)) + push!(loss_sq, dot(Δ_proj, Δ_proj)) + end + mesh, loss_lin, loss_sq +end + +function circles_in_triangle_shape() + # initialize the partial gram matrix for a sphere inscribed in a regular + # tetrahedron + J = Int64[] + K = Int64[] + values = BigFloat[] + for j in 1:8 + for k in 1:8 + filled = false + if j == k + push!(values, 1) + filled = true + elseif (j == 1 || k == 1) + push!(values, 0) + filled = true + elseif (j == 2 || k == 2) + push!(values, -1) + filled = true + end + #=elseif (j <= 5 && j != 2 && k == 9 || k == 9 && k <= 5 && k != 2) + push!(values, 0) + filled = true + end=# + if filled + push!(J, j) + push!(K, k) + end + end + end + append!(J, [6, 4, 6, 5, 7, 5, 7, 3, 8, 3, 8, 4]) + append!(K, [4, 6, 5, 6, 5, 7, 3, 7, 3, 8, 4, 8]) + append!(values, fill(-1, 12)) + + # plot loss along a slice + loss_lin = [] + loss_sq = [] + mesh = range(0.99, 1.01, 101) + for t in mesh + L = hcat( + Engine.plane(BigFloat[0, 0, 1], BigFloat(0)), + Engine.sphere(BigFloat[0, 0, 0], BigFloat(t)), + Engine.plane(BigFloat[1, 0, 0], BigFloat(1)), + Engine.plane(BigFloat[cos(2pi/3), sin(2pi/3), 0], BigFloat(1)), + Engine.plane(BigFloat[cos(-2pi/3), sin(-2pi/3), 0], BigFloat(1)), + Engine.sphere(4//3*BigFloat[-1, 0, 0], BigFloat(1//3)), + Engine.sphere(4//3*BigFloat[cos(-pi/3), sin(-pi/3), 0], BigFloat(1//3)), + Engine.sphere(4//3*BigFloat[cos(pi/3), sin(pi/3), 0], BigFloat(1//3)) + ) + Δ_proj = Engine.proj_diff(gram, L'*Engine.Q*L) + push!(loss_lin, norm(Δ_proj)) + push!(loss_sq, dot(Δ_proj, Δ_proj)) + end + mesh, loss_lin, loss_sq +end \ No newline at end of file diff --git a/engine-proto/gram-test/circles-in-triangle.jl b/engine-proto/gram-test/circles-in-triangle.jl index 2bca4f2..b031711 100644 --- a/engine-proto/gram-test/circles-in-triangle.jl +++ b/engine-proto/gram-test/circles-in-triangle.jl @@ -1,20 +1,28 @@ include("Engine.jl") using SparseArrays -using AbstractAlgebra -using PolynomialRoots # initialize the partial gram matrix for a sphere inscribed in a regular # tetrahedron J = Int64[] K = Int64[] values = BigFloat[] -for j in 1:8 - for k in 1:8 +for j in 1:9 + for k in 1:9 filled = false if j == k - push!(values, 1) + push!(values, j < 9 ? 1 : 0) filled = true + elseif (j == 9) + if (k <= 5 && k != 2) + push!(values, 0) + filled = true + end + elseif (k == 9) + if (j <= 5 && j != 2) + push!(values, 0) + filled = true + end elseif (j == 1 || k == 1) push!(values, 0) filled = true @@ -47,7 +55,6 @@ gram = sparse(J, K, values) ## guess = Engine.rand_on_shell(fill(BigFloat(-1), 8)) # set initial guess -#= guess = hcat( Engine.plane(BigFloat[0, 0, 1], BigFloat(0)), Engine.sphere(BigFloat[0, 0, 0], BigFloat(1//2)), @@ -56,9 +63,10 @@ guess = hcat( Engine.plane(BigFloat[cos(-2pi/3), sin(-2pi/3), 0], BigFloat(1)), Engine.sphere(BigFloat[-1, 0, 0], BigFloat(1//5)), Engine.sphere(BigFloat[cos(-pi/3), sin(-pi/3), 0], BigFloat(1//5)), - Engine.sphere(BigFloat[cos(pi/3), sin(pi/3), 0], BigFloat(1//5)) + Engine.sphere(BigFloat[cos(pi/3), sin(pi/3), 0], BigFloat(1//5)), + BigFloat[0, 0, 0, 1, 1] ) -=# +#= guess = hcat( Engine.plane(BigFloat[0, 0, 1], BigFloat(0)), Engine.sphere(BigFloat[0, 0, 0], BigFloat(0.9)), @@ -67,13 +75,22 @@ guess = hcat( Engine.plane(BigFloat[cos(-2pi/3), sin(-2pi/3), 0], BigFloat(1)), Engine.sphere(4//3*BigFloat[-1, 0, 0], BigFloat(1//3)), Engine.sphere(4//3*BigFloat[cos(-pi/3), sin(-pi/3), 0], BigFloat(1//3)), - Engine.sphere(4//3*BigFloat[cos(pi/3), sin(pi/3), 0], BigFloat(1//3)) + Engine.sphere(4//3*BigFloat[cos(pi/3), sin(pi/3), 0], BigFloat(1//3)), + BigFloat[0, 0, 0, 1, 1] ) +=# -# complete the gram matrix using gradient descent -L, history = Engine.realize_gram(gram, guess, max_descent_steps = 200) +# 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) completed_gram = L'*Engine.Q*L println("Completed Gram matrix:\n") display(completed_gram) -println("\nSteps: ", size(history.stepsize, 1)) -println("Loss: ", history.scaled_loss[end], "\n") \ No newline at end of file +println( + "\nSteps: ", + size(history.scaled_loss, 1), + " + ", size(history_pol.scaled_loss, 1), + " + ", size(history_pol2.scaled_loss, 1) +) +println("Loss: ", history_pol2.scaled_loss[end], "\n") \ No newline at end of file diff --git a/engine-proto/gram-test/overlapping-pyramids.jl b/engine-proto/gram-test/overlapping-pyramids.jl index 97fc119..51606e1 100644 --- a/engine-proto/gram-test/overlapping-pyramids.jl +++ b/engine-proto/gram-test/overlapping-pyramids.jl @@ -47,13 +47,14 @@ guess = hcat( Engine.rand_on_shell(fill(BigFloat(-1), 2)) ) -# complete the gram matrix using gradient descent -L, history = Engine.realize_gram(gram, guess) +# 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) completed_gram = L'*Engine.Q*L println("Completed Gram matrix:\n") display(completed_gram) -println("\nSteps: ", size(history.stepsize, 1)) -println("Loss: ", history.scaled_loss[end], "\n") +println("\nSteps: ", size(history.scaled_loss, 1), " + ", size(history_pol.scaled_loss, 1)) +println("Loss: ", history_pol.scaled_loss[end], "\n") # === algebraic check === diff --git a/engine-proto/gram-test/sphere-in-tetrahedron.jl b/engine-proto/gram-test/sphere-in-tetrahedron.jl index 3730390..dc5588c 100644 --- a/engine-proto/gram-test/sphere-in-tetrahedron.jl +++ b/engine-proto/gram-test/sphere-in-tetrahedron.jl @@ -1,8 +1,6 @@ include("Engine.jl") using SparseArrays -using AbstractAlgebra -using PolynomialRoots using Random # initialize the partial gram matrix for a sphere inscribed in a regular @@ -35,10 +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 using gradient descent -L, history = Engine.realize_gram(gram, guess) +# complete the gram matrix using Newton's method +L, history = Engine.realize_gram_newton(gram, guess) completed_gram = L'*Engine.Q*L println("Completed Gram matrix:\n") display(completed_gram) -println("\nSteps: ", size(history.stepsize, 1)) +println("\nSteps: ", size(history.scaled_loss, 1)) println("Loss: ", history.scaled_loss[end], "\n") \ No newline at end of file