I am trying to figure out how generally to go about solving an energy minimization problem such as the one below. I am just not used to dealing with large systems, as most of the maths that I have learned has been much more "small-scale".
This is part of implementing a solver for mass-spring systems (outlined in this paper) that involves alternating between solving for the spring directions, and then the node positions. This technique is referred to as "alternating optimization", or "Block Coordinate Descent" by the authors.
Initially the node positions are estimated, and the directions are calculated, following which the directions are fixed and the positions of each node are computed.
The paper fairly succinctly boils the problem down to the following equation: $$ \frac{1}{2}x^T(M+h^2L)x-h^2x^TJd+x^Tb $$
Where:
- the system consists of $m$ nodes and $s$ springs
- $h$ is the timestep(constant)
- $x$ is a column vector of length $3m$ (or $3m*1$ matrix) representing the node positions,
- $(M+h^2L)$ is a matrix of size $3m*3m$ (which, in this case, is symmetric, positive semi-definite, and does not change after being initialized.),
- $J$ is a matrix of size $3m*3s$, and
- $d$ is a column vector of size $3s$ representing the spring directions. (this particular equation is equation 14 in section 4 of the paper)
- $b$ is a column vector of length $3m$, which represents the xyz components of the external forces acting on each node
I am at a stage where I can generate all of vectors and matrices in code based on some arbitrary mesh consisting of nodes and links, but I am having trouble figuring out how I would go about solving for $d$ and $x$. I have an initial guess for $x$ (as outlined towards the ends of both section 3 and 4), but I don't have much experience solving entire systems like this.
My main question has two parts:
- is there some generic procedure one would use to go about solving such a system?
- if not, could some resources on the topic be suggested? Even some search terms more specific than "energy minimisation" would be of great use!
Ultimately I just don't know where to start, and the paper seems to simply say "starting with an initial guess for $x$, first we compute the optimal $d$", but I don't know what "compute the optimal $d$" entails.
Similarly, I wouldn't know where to start with computing the optimal values for x. I assume this would involve differentiating w.r.t. x and finding the values when the derivative equals zero but again, I have no idea how to apply this to entire matrices of values, as opposed to single values.
I am aware that each term in the above equation should evaluate to a scalar value, but I don't know how I would be able to obtain a value for each row/cell in the column vectors $x$ and $d$.
NOTE: I believe there are one or two minor errors in the paper, such as $d ∈ R^{2s}$ which, I believe, should be $d ∈ R^{3s}$, since we are dealing with springs in 3 dimensions.
According to your comment above you're looking for the following:
\begin{align} \arg \min_{x} \left\{ \frac{1}{2} {x}^{T} \left( M + {h}^{2} L \right) x - {h}^{2} {x}^{T} J d + {x}^{T} b \right\} \end{align}
Rewriting the above yields:
\begin{align} \arg \min_{x} \left\{ \frac{1}{2} {x}^{T} \left( M + {h}^{2} L \right) x - {\left( {h}^{2} J d + b \right)}^{T} x \right\} \end{align}
Now, according to your description $ M + {h}^{2} L $ is PSD Matrix hence the above is a simple quadratic form of a convex problem:
\begin{align} \arg \min_{x} \left\{ \frac{1}{2} {x}^{T} \left( M + {h}^{2} L \right) x - {\left( {h}^{2} J d + b \right)}^{T} x \right\} & = \arg \min_{x} \left\{ \frac{1}{2} {x}^{T} A x - {c}^{T} x \right\} \\ & = {A}^{-1} c \\ & = {\left( M + {h}^{2} L \right)}^{-1} \left( {h}^{2} J d + b \right) \end{align}