I don't know much about tracking implicit surfaces, so I'm just going to start with the optimization problem and go from there.
The optimization problem is, at the core, nonlinear least squares, and can therefore be solved efficiently by the Gauss-Newton method, wherein at each step one linearizes the problem at the current guess, solves the linear least squares problem for this linearization, then uses the solution as the initial guess for the next step.
At each step, the linear least squares problem can be recast as solving one large sparse linear system. Figuring out what the coefficient matrix is for this linear system is a little mindbending, but doable with some matrix vectorization tricks (explanation below).
Once the matrix is constructed, the linear system can be solved efficiently with any good sparse-direct factorization technique up to quite large problem sizes. Furthermore, this matrix is physically analogous to the matrices that arise when trying to find the minimum energy configuration of a spring network, and so if the problem is extremely huge (too big for a sparse-direct factorization to fit in memory) it should still be amenable to multigrid.
This optimization problem has the interesting property that the data misfit term is linear (and therefore stays the same regardless of the iteration), whereas the regularization term is nonlinear (and therefore changes from iteration to iteration). Usually it is exactly the opposite, where the data misfit is nonlinear but the regularization is linear.
Optimization problem
In what follows, let $N$ be the total number of grid vertices, and let $M$ be the total number of target points $p'$.
If I read everything correctly, the optimization problem can be recast as:
$$\min_{V'} \frac{1}{2}||\hat{Q}'(V') - P'||^2_\text{fro} + \frac{1}{2}||D~m(V'-V)||^2_\text{fro}$$
In this formula,
- $||\cdot||_\text{fro}$ is the Frobenius matrix norm (square root of sum of squares of all the entries of the matrix),
- $V$ is the $3$-by-$N$ matrix where each column is one of the vertices $v$ in the initial configuration.
- $V'$ is the $3$-by-$N$ matrix where each column is one of the vertices $v'$ of the moved configuration.
- $P'$ is the $3$-by-$M$ matrix where each column is one of the target points $p'$.
- $\hat{Q}'(V')$ is the $3$-by-$M$ matrix of approximate target points, $q' \approx p'$, based on a given choice of moved vertices $V'$.
- $D$ is the $M$-by-$N$ graph finite difference matrix that takes in a vector of values on the nodes of a graph, and outputs a list of differences across the edges. The $j$'th entry of the output is the difference between the value at the head vertex and the value at the tail vertex for edge $e_j$. I.e., the $i$'th row of this matrix has $1$ in the column corresponding to the "head" vertex for edge $e_i$, $-1$ in the column corresponding to the "tail" vertex for edge $e_i$, and zeros elsewhere.
- $m:\mathbb{R}^{3\times N}\rightarrow \mathbb{R^N}$ is the columnwise length measuring function. I.e, if $X$ is a matrix with columns $x_i$, then the $i$'th element of $m(X)$ is $||x_i||$.
Data misfit term
Suppose the $j$'th approximate vertex, $q_j$, lies on an edge, where the "tail" of the edge is a moved vertex $v'_\text{tail}$, and the "head" of the edge is a moved vertex $v'_\text{head}$. Then from the problem statement,
$$q = (1-\alpha_j)v'_\text{head} + \alpha_j v'_\text{tail}.$$
Then, applying this transformation to all vertices at once, we have:
$$Q'(V') = (1 - A).*(V' E_\text{heads}) + A.*(V' E_\text{tails}),$$
where
- "$.*$" is MATLAB-style notation for the elementwise (Hadamard) matrix product.
- $A$ is the matrix where each column is constant, and the $j$'th column contains copies of $\alpha_j$ repeated vertically,
- $1$ represents the matrix of all ones, of the same size as $A$ (not the identity matrix),
- $E_\text{heads}$ and $E_\text{tails}$ are the column selection matrices, that, when acting on the matrix $V'$ from the right, generate the matrices whose $j$'th columns are the "head" and "tail" vectors for $q_j$, respectively. More details on these column selection matrices later..
Element-wise multiplication of matrices is the same as multipying the vectorized form of one of the matrices by a big diagonal matrix whose diagonal entries come from the other matrix. Thus if we vectorizes (unrolls) the matrix $Q'(V')$ into one long vector, "$\text{vec}(Q'(V')$", we can write it as follows:
$$\text{vec}(Q'(V')) = (I_{3M}-D)~\text{vec}(V' E_\text{heads}) + D~\text{vec}(V' E_\text{tails}),$$
where $D:=\text{diag}(\text{vec}(A))$ is a very big diagonal matrix whose diagonal entries are the entries of $A$, and the notation $I_{k}$ denotes a $k$-by-$k$ identity matrix.
Furthermore, recall the following relationship between vectorization $\text{vec}(\cdot)$ and Kronecker products $\otimes$:
$$\text{vec}(ABC) = (C^T \otimes A)\text{vec}(B).$$
Applying this identity to our expression above, we have:
$$\text{vec}(Q'(V')) = (I_{3M}-D)~(E_\text{heads}^T \otimes I_3)~\text{vec}(V') + D~(E_\text{tails}^T \otimes I_3)~\text{vec}(V'),$$
or
$$\text{vec}(Q'(V')) = J~ \text{vec}(V'),$$
where
$$\boxed{J:=(I_{3M}-D)~(E_\text{heads}^T \otimes I_3) + D~(E_\text{tails}^T \otimes I_3)}$$
is the matrix that describes the vectorized transformation from all moved $v'$ vertices to all approximate target points $q'$. This matrix is extremely large, but also extremely sparse. (Also very easy to construct in MATLAB from built-in functions)
Regularization term
Recall that the regularization term is:
$$\frac{1}{2}||D~m(V'-V)||^2_\text{fro},$$
where, $m(\cdot)$ is the function that takes in a $3$-by-$N$ matrix, and outputs the vector containing the lengths of each column of the input matrix. I.e.,
$$m(X) = \begin{bmatrix}
||x_1|| \\
||x_2|| \\
\vdots \\
||x_N|| \\
\end{bmatrix}$$
We seek to linearize this term about some initial guess, $V'_0$, for the moved point $V'$.
Using the formula for the derivative of the norm, we can compute the the derivative of $m$ at point $Y$ in direction $\delta X$ as follows:
\begin{align}
\left.\frac{dm}{dX}\right|_{Y}~\delta X &= \begin{bmatrix}
\frac{y_1^T \delta x_1}{||y_1||} \\
\frac{y_2^T \delta x_2}{||y_2||} \\
\vdots \\
\frac{y_N^T \delta x_N}{||y_N||} \\
\end{bmatrix} \\
&= \begin{bmatrix}
\frac{y_1^T}{||y_1||} \\
& \frac{y_2^T}{||y_2||} \\
& & \ddots \\
& & & \frac{y_N^T}{||y_N||}
\end{bmatrix}
\begin{bmatrix}
\delta x_1 \\
\delta x_2 \\
\vdots \\
\delta x_N
\end{bmatrix} \\
&= M(Y)~\text{vec}(\delta X),
\end{align}
where $M(Y)$ is the big $N$-by-$3N$ matrix shown in the middle formula. I.e.,
$$M(Y)=\text{blockdiag}\left(\frac{y_1^T}{||y_1||}, \frac{y_2^T}{||y_2||}, \dots, \frac{y_N^T}{||y_N||}\right).$$
Hence, the regularization function can be linearized about an initial point $V'_0 \approx V'$ as follows,
$$D~m(V'-V) = D~m(V_0' - V) + D~M(V_0' - V)~\text{vec}(V' - V'_0) + O(||V' - V'_0||^2),$$
Defining
$$\boxed{R(V, V_0') := D~M(V_0' - V)}$$
and
$$\boxed{r_0(V, V'_0) := D~m(V_0' - V) - R(V, V_0')~\text{vec}(V'_0),}$$
the linearization of the regularization term can be written as follows:
$$D~m(V'-V) = r_0(V, V_0') + R(V, V_0')~\text{vec}(V') + O(||V' - V'_0||^2).$$
Linearized optimization problem
Hence, the original optimization problem, when the regularization is linearized about a point $V'_0$, can be recast in vectorized form as follows:
$$\min_{\text{vec}(V')} \frac{1}{2}||J~\text{vec}(V') - \text{vec}(P')||^2 + \frac{1}{2}||R~\text{vec}(V') + r_0||^2,$$
where we abbreviate $R:=R(V, V_0')$, and $r_0:=r_0(V, V_0')$. This linear least squares problem can be solved by solving the normal equations system:
$$\boxed{(J^TJ + R^TR)~\text{vec}(V') = J^T~\text{vec}(P') - R^T~r_0.}$$
So, you can form solve this sparse system for $\text{vec}(V')$, then reshape the result from a vector into a matrix, and then extract the columns of this matrix as the desired vertices $v'$.
This will give you the optimal moved vertices based on the linearization of the regularization about point $V'_0$. To solve the problem with general nonlinear regularization, you can iterate: linearize about an initial guess, solve the linearized problem, use the solution as the next point to linearize about, solve that problem, and repeat until convergence. This is the Gauss-Newton method.
A side note of interest is that the graph Laplacian $D^TD$ appears within the core of the $R^TR$, so this regularization could be interpreted as some form of grouped discrete smoothing regularization.
Heads and tails selection matrices
Above we glossed over the construction of the "heads" and "tails" selection matrices $E_\text{heads}$ and $E_\text{tails}$, so more details are provided here.
Let $e_\text{heads}$ and $e_\text{tails}$ be length-$M$ vectors representing the "head-for-$p$" map and "tail-for-$p$" maps. That is, if $p_j$ lies on the edge going from vertex $v_a$ to vertex $v_b$, then the $j$'th entry of $e_\text{heads}$ is the index $b$, and the $j$'th entry of $e_\text{tails}$ is the index $a$.
For example, suppose the first point lies on the edge going from vertex $v_3$ to $v_8$, the second point lies on the edge going from vertex $v_5$ to $v_6$, and the third point lies on the edge going from $v_9$ to $v_5$ (this is a made up example to illustrate the point). Schematically,
\begin{align}
v_3 &\overset{p_1}{\rightarrow} v_8 \\
v_5 &\overset{p_2}{\rightarrow} v_6 \\
v_9 &\overset{p_3}{\rightarrow} v_5.
\end{align}
Then the "edge-to-head" and "edge-to-tail" vectors would start:
\begin{align}
e_\text{tails}=&\begin{bmatrix}8, & 6, & 5, & \dots\end{bmatrix} \\
e_\text{heads}=&\begin{bmatrix}3, & 5, & 9, & \dots\end{bmatrix}.
\end{align}
We want to find the selection matrices $E_\text{heads}$ and $E_\text{tails}$ such that:
$$\begin{bmatrix}v_1 & v_2 & \dots & v_N\end{bmatrix} E_\text{heads}
=
\begin{bmatrix}v_{e_\text{heads}(1)} & v_{e_\text{heads}(2)} & \dots & v_{e_\text{heads}(M)}\end{bmatrix},$$
and
$$\begin{bmatrix}v_1 & v_2 & \dots & v_N\end{bmatrix} E_\text{tails}
=
\begin{bmatrix}v_{e_\text{tails}(1)} & v_{e_\text{tails}(2)} & \dots & v_{e_\text{tails}(M)}\end{bmatrix}.$$
The matrix $E_\text{heads}$ that does this is the sparse $N$-by-$M$ matrix, where the $j$'th column has a $1$ in row $e_\text{heads}(j)$ (and likewise for $E_\text{tails}$, except with the row indices drawn from $e_\text{tails}$).
For example, consider the following product of a matrix whose columns are the vectors $x_1$, $x_2$, and $x_3$, with a selection matrix:
$$\begin{bmatrix}
x_1 & x_2 & x_3
\end{bmatrix}
\begin{bmatrix}
& & 1\\
& & & 1 \\
1 & 1 & &
\end{bmatrix}
=
\begin{bmatrix}
x_3 & x_3 & x_1 & x_2
\end{bmatrix}.$$
