From 2c55a63a6ff959bb389af93e5677815cf4b87117 Mon Sep 17 00:00:00 2001 From: Aaron Fenyes Date: Fri, 6 Dec 2024 14:35:30 -0800 Subject: [PATCH] Engine: Find the tangent space of the solution variety At the end of the realization routine, use the computed Hessian to find the tangent space of the solution variety, and return it alongside the realization. Since altering the constraints can change the tangent space without changing the solution, we compute the tangent space even when the guess passed to the realization routine is already a solution. --- app-proto/examples/irisawa-hexlet.rs | 2 +- app-proto/examples/point-on-sphere.rs | 2 +- app-proto/examples/three-spheres.rs | 2 +- app-proto/src/assembly.rs | 12 ++- app-proto/src/engine.rs | 131 +++++++++++++++++++++++--- 5 files changed, 129 insertions(+), 20 deletions(-) diff --git a/app-proto/examples/irisawa-hexlet.rs b/app-proto/examples/irisawa-hexlet.rs index fc14f91..2bc94c0 100644 --- a/app-proto/examples/irisawa-hexlet.rs +++ b/app-proto/examples/irisawa-hexlet.rs @@ -2,7 +2,7 @@ use dyna3::engine::{Q, irisawa::realize_irisawa_hexlet}; fn main() { const SCALED_TOL: f64 = 1.0e-12; - let (config, success, history) = realize_irisawa_hexlet(SCALED_TOL); + let (config, _, success, history) = realize_irisawa_hexlet(SCALED_TOL); print!("\nCompleted Gram matrix:{}", config.tr_mul(&*Q) * &config); if success { println!("Target accuracy achieved!"); diff --git a/app-proto/examples/point-on-sphere.rs b/app-proto/examples/point-on-sphere.rs index 0e4765a..13040e5 100644 --- a/app-proto/examples/point-on-sphere.rs +++ b/app-proto/examples/point-on-sphere.rs @@ -18,7 +18,7 @@ fn main() { ]); let frozen = [(3, 0)]; println!(); - let (config, success, history) = realize_gram( + let (config, _, success, history) = realize_gram( &gram, guess, &frozen, 1.0e-12, 0.5, 0.9, 1.1, 200, 110 ); diff --git a/app-proto/examples/three-spheres.rs b/app-proto/examples/three-spheres.rs index d348b18..19acfd1 100644 --- a/app-proto/examples/three-spheres.rs +++ b/app-proto/examples/three-spheres.rs @@ -21,7 +21,7 @@ fn main() { ]) }; println!(); - let (config, success, history) = realize_gram( + let (config, _, success, history) = realize_gram( &gram, guess, &[], 1.0e-12, 0.5, 0.9, 1.1, 200, 110 ); diff --git a/app-proto/src/assembly.rs b/app-proto/src/assembly.rs index fb5bbf7..0d63e45 100644 --- a/app-proto/src/assembly.rs +++ b/app-proto/src/assembly.rs @@ -5,7 +5,7 @@ use std::{collections::BTreeSet, sync::atomic::{AtomicU64, Ordering}}; use sycamore::prelude::*; use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */ -use crate::engine::{realize_gram, PartialMatrix}; +use crate::engine::{realize_gram, ConfigSubspace, PartialMatrix}; // the types of the keys we use to access an assembly's elements and constraints pub type ElementKey = usize; @@ -109,7 +109,6 @@ impl Element { } } } - #[derive(Clone)] pub struct Constraint { @@ -127,6 +126,9 @@ pub struct Assembly { pub elements: Signal>, pub constraints: Signal>, + // solution variety tangent space + pub tangent: Signal, + // indexing pub elements_by_id: Signal> } @@ -136,6 +138,7 @@ impl Assembly { Assembly { elements: create_signal(Slab::new()), constraints: create_signal(Slab::new()), + tangent: create_signal(ConfigSubspace::zero()), elements_by_id: create_signal(FxHashMap::default()) } } @@ -247,7 +250,7 @@ impl Assembly { } // look for a configuration with the given Gram matrix - let (config, success, history) = realize_gram( + let (config, tangent, success, history) = realize_gram( &gram, guess, &[], 1.0e-12, 0.5, 0.9, 1.1, 200, 110 ); @@ -271,6 +274,9 @@ impl Assembly { |rep| rep.set_column(0, &config.column(elt.column_index)) ); } + + // save the tangent space + self.tangent.set_silent(tangent); } } } \ No newline at end of file diff --git a/app-proto/src/engine.rs b/app-proto/src/engine.rs index 285a9c6..6a62579 100644 --- a/app-proto/src/engine.rs +++ b/app-proto/src/engine.rs @@ -1,5 +1,5 @@ use lazy_static::lazy_static; -use nalgebra::{Const, DMatrix, DVector, Dyn}; +use nalgebra::{Const, DMatrix, DVector, DVectorView, Dyn, SymmetricEigen}; use web_sys::{console, wasm_bindgen::JsValue}; /* DEBUG */ // --- elements --- @@ -85,6 +85,47 @@ impl PartialMatrix { } } +// --- configuration subspaces --- + +pub struct ConfigSubspace(Vec>); + +impl ConfigSubspace { + pub fn zero() -> ConfigSubspace { + ConfigSubspace(Vec::new()) + } + + // approximate the kernel of a symmetric endomorphism of the configuration + // space for `assembly_dim` elements. we consider an eigenvector to be part + // of the kernel if its eigenvalue is smaller than the constant `THRESHOLD` + fn symmetric_kernel(a: DMatrix, assembly_dim: usize) -> ConfigSubspace { + const ELEMENT_DIM: usize = 5; + const THRESHOLD: f64 = 1.0e-9; + let eig = SymmetricEigen::new(a); + let eig_vecs = eig.eigenvectors.column_iter(); + let eig_pairs = eig.eigenvalues.iter().zip(eig_vecs); + let basis = eig_pairs.filter_map( + |(λ, v)| (λ.abs() < THRESHOLD).then_some( + Into::>::into( + v.reshape_generic(Dyn(ELEMENT_DIM), Dyn(assembly_dim)) + ) + ) + ); + ConfigSubspace(basis.collect()) + } + + // find the projection onto this subspace of the motion where the element + // with the given column index has velocity `v` + /* TO DO */ + // for the zero subspace, this method's behavior doesn't match its name: it + // panics rather than returning zero + fn proj(&self, v: &DVectorView, column_index: usize) -> DMatrix { + let ConfigSubspace(basis) = self; + basis.into_iter().map( + |b| b.column(column_index).dot(&v) * b + ).sum() + } +} + // --- descent history --- pub struct DescentHistory { @@ -181,7 +222,7 @@ pub fn realize_gram( reg_scale: f64, max_descent_steps: i32, max_backoff_steps: i32 -) -> (DMatrix, bool, DescentHistory) { +) -> (DMatrix, ConfigSubspace, bool, DescentHistory) { // start the descent history let mut history = DescentHistory::new(); @@ -201,12 +242,8 @@ pub fn realize_gram( // use Newton's method with backtracking and gradient descent backup let mut state = SearchState::from_config(gram, guess); + let mut hess = DMatrix::zeros(element_dim, assembly_dim); for _ in 0..max_descent_steps { - // stop if the loss is tolerably low - history.config.push(state.config.clone()); - history.scaled_loss.push(state.loss / scale_adjustment); - if state.loss < tol { break; } - // find the negative gradient of the loss function let neg_grad = 4.0 * &*Q * &state.config * &state.err_proj; let mut neg_grad_stacked = neg_grad.clone().reshape_generic(Dyn(total_dim), Const::<1>); @@ -229,7 +266,7 @@ pub fn realize_gram( hess_cols.push(deriv_grad.reshape_generic(Dyn(total_dim), Const::<1>)); } } - let mut hess = DMatrix::from_columns(hess_cols.as_slice()); + hess = DMatrix::from_columns(hess_cols.as_slice()); // regularize the Hessian let min_eigval = hess.symmetric_eigenvalues().min(); @@ -249,6 +286,11 @@ pub fn realize_gram( hess[(k, k)] = 1.0; } + // stop if the loss is tolerably low + history.config.push(state.config.clone()); + history.scaled_loss.push(state.loss / scale_adjustment); + if state.loss < tol { break; } + // compute the Newton step /* we need to either handle or eliminate the case where the minimum @@ -256,7 +298,7 @@ pub fn realize_gram( singular. right now, this causes the Cholesky decomposition to return `None`, leading to a panic when we unrap */ - let base_step_stacked = hess.cholesky().unwrap().solve(&neg_grad_stacked); + let base_step_stacked = hess.clone().cholesky().unwrap().solve(&neg_grad_stacked); let base_step = base_step_stacked.reshape_generic(Dyn(element_dim), Dyn(assembly_dim)); history.base_step.push(base_step.clone()); @@ -269,10 +311,16 @@ pub fn realize_gram( state = better_state; history.backoff_steps.push(backoff_steps); }, - None => return (state.config, false, history) + None => return (state.config, ConfigSubspace::zero(), false, history) }; } - (state.config, state.loss < tol, history) + let success = state.loss < tol; + let tangent = if success { + ConfigSubspace::symmetric_kernel(hess, assembly_dim) + } else { + ConfigSubspace::zero() + }; + (state.config, tangent, success, history) } // --- tests --- @@ -291,7 +339,7 @@ pub mod irisawa { use super::*; - pub fn realize_irisawa_hexlet(scaled_tol: f64) -> (DMatrix, bool, DescentHistory) { + pub fn realize_irisawa_hexlet(scaled_tol: f64) -> (DMatrix, ConfigSubspace, bool, DescentHistory) { let gram = { let mut gram_to_be = PartialMatrix::new(); for s in 0..9 { @@ -399,7 +447,7 @@ mod tests { fn irisawa_hexlet_test() { // solve Irisawa's problem const SCALED_TOL: f64 = 1.0e-12; - let (config, _, _) = realize_irisawa_hexlet(SCALED_TOL); + let (config, _, _, _) = realize_irisawa_hexlet(SCALED_TOL); // check against Irisawa's solution let entry_tol = SCALED_TOL.sqrt(); @@ -409,6 +457,61 @@ mod tests { } } + #[test] + fn tangent_test() { + const SCALED_TOL: f64 = 1.0e-12; + const ELEMENT_DIM: usize = 5; + const ASSEMBLY_DIM: usize = 3; + let gram = { + let mut gram_to_be = PartialMatrix::new(); + for j in 0..3 { + for k in j..3 { + gram_to_be.push_sym(j, k, if j == k { 1.0 } else { -1.0 }); + } + } + gram_to_be + }; + let guess = DMatrix::from_columns(&[ + sphere(0.0, 0.0, 0.0, -2.0), + sphere(0.0, 0.0, 1.0, 1.0), + sphere(0.0, 0.0, -1.0, 1.0) + ]); + let frozen: [_; 5] = std::array::from_fn(|k| (k, 0)); + let (config, tangent, success, history) = realize_gram( + &gram, guess.clone(), &frozen, + SCALED_TOL, 0.5, 0.9, 1.1, 200, 110 + ); + assert_eq!(config, guess); + assert_eq!(success, true); + assert_eq!(history.scaled_loss.len(), 1); + + // confirm that the tangent space has dimension five or less + let ConfigSubspace(ref tangent_basis) = tangent; + assert_eq!(tangent_basis.len(), 5); + + // confirm that the tangent space contains all the motions we expect it + // to. since we've already bounded the dimension of the tangent space, + // this confirms that the tangent space is what we expect it to be + let tangent_motions = vec![ + basis_matrix((0, 1), ELEMENT_DIM, ASSEMBLY_DIM), + basis_matrix((1, 1), ELEMENT_DIM, ASSEMBLY_DIM), + basis_matrix((0, 2), ELEMENT_DIM, ASSEMBLY_DIM), + basis_matrix((1, 2), ELEMENT_DIM, ASSEMBLY_DIM), + DMatrix::::from_column_slice(ELEMENT_DIM, 3, &[ + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, -1.0, -0.25, -1.0, + 0.0, 0.0, -1.0, 0.25, 1.0 + ]) + ]; + let tol_sq = ((ELEMENT_DIM * ASSEMBLY_DIM) as f64) * SCALED_TOL * SCALED_TOL; + for motion in tangent_motions { + let motion_proj: DMatrix<_> = motion.column_iter().enumerate().map( + |(k, v)| tangent.proj(&v, k) + ).sum(); + assert!((motion - motion_proj).norm_squared() < tol_sq); + } + } + // at the frozen indices, the optimization steps should have exact zeros, // and the realized configuration should match the initial guess #[test] @@ -428,7 +531,7 @@ mod tests { ]); let frozen = [(3, 0), (3, 1)]; println!(); - let (config, success, history) = realize_gram( + let (config, _, success, history) = realize_gram( &gram, guess.clone(), &frozen, 1.0e-12, 0.5, 0.9, 1.1, 200, 110 );