The following problem was posted on the Julia Discourse:
So I have three magnetic resonance (MR) images: $X, Y, Z$. They come from three different pulse sequences which interact with NMR relaxation parameters in known ways. That is, we have equations which given the underlying NMR parameters (specifically, T1, T2 and PD) and the scanner parameters $p^X, p^Y, p^Z$ corresponding to each image $X,Y,Z$, we can model the output intensity value observed in the acquired image at some voxel [$(i, j, k)$, e.g., $X_{i,j,k}$, since we are in 3D].
However, we don’t know the $T_1$, $T_2$ and $PD$ parameters but we do have $X,Y,Z$ and we can separately solve for each $p^{(\cdot)}$. If you haven’t followed any of the previous mediocre explanation, then the problem simply reduces to the following system of equations, for which I am trying to solve $T_1$, $T_2$, and $P_D$ where everything else is given.
$$
\begin{align} \log(X_{i,j,k}) &= \log(PD) + p_0^X + p_1^X T_1 + p_2^X T_1^2 \\ \log(Y_{i,j,k}) &= \log(PD) + p_0^Y + p_1^Y T_1 + p_2^Y / T_2 \\ \log(Z_{i,j,k}) &= \log(PD) + p_0^Z + p_1^Z T_1 + p_2^Z / T_2 \\ \end{align} $$
In previous iterations of this project, the person who built the code used non-linear least squares (in MATLAB, FWIW) to solve for the values of $T_1$, $T_2$ and $PD$. However, the person who built the code left our lab, and with that has gone the knowledge of why that was chosen.
Previously I was using python’s scipy implementation of least squares to solve the problem at every voxel; however, even with parallelization this takes hours.
Since these are polynomial equations we can use homotopy continuation to solve the problem. Let's setup the system first:
using HomotopyContinuation, DynamicPolynomials
@polyvar x y z PD p[1:3,1:3] T[1:2]
f1 = PD + p[1,1] + p[1,2]*T[1] + p[1,3] * T[1]^2 - x
# note that we eliminate the denominator in equation 2 and 3
f2 = (PD + p[2,1] + p[2,2]*T[1])*T[2] + p[2,3] - y * T[2]
f3 = (PD + p[3,1] + p[3,2]*T[1])*T[2] + p[3,3] - z * T[2]
f = [f1, f2, f3]
3-element Array{Polynomial{true,Int64},1}: p₁₋₃T₁² + p₁₋₂T₁ - x + PD + p₁₋₁ p₂₋₂T₁T₂ - yT₂ + PDT₂ + p₂₋₁T₂ + p₂₋₃ p₃₋₂T₁T₂ - zT₂ + PDT₂ + p₃₋₁T₂ + p₃₋₃
Since we need to solve the same instance many many times it makes sense to first compute the generic number of solutions of this polynomial system and then to use a parameter homotopy.
We need some generic points to generate a generic instance
params = [vec(p); x; y; z]
generic_params = randn(ComplexF64, length(params)) # Note we should sample them randomly
f_generic = subs.(f, Ref(params => generic_params))
3-element Array{Polynomial{true,Complex{Float64}},1}: (0.06560018192766423 + 1.226509609243317im)T₁² + PD + (-0.40076878384840564 + 0.3936976384006835im)T₁ + (-0.4034430867072931 + 2.130445251499418im) PDT₂ + (0.17548682027400092 + 1.0868781469710116im)T₁T₂ + (1.0068779513453254 - 0.9349141017423815im)T₂ + (0.6579055458504889 + 0.11880590069556508im) PDT₂ + (-0.6911103187839271 - 1.1844653755707422im)T₁T₂ + (0.216353296242008 - 0.3087205812565972im)T₂ + (-1.3312172733617187 + 1.9140562766862028im)
generic_result = solve(f_generic)
AffineResult with 8 tracked paths ================================== • 2 non-singular finite solutions (0 real) • 0 singular finite solutions (0 real) • 6 solutions at infinity • 0 failed paths • random seed: 808747
We see that the system has only 2 solutions (instead the maximal possible 8 solutions).
generic_solutions = solutions(generic_result)
2-element Array{Array{Complex{Float64},1},1}: [0.100221+0.505434im, -0.256837+1.29393im, 0.809632-0.787296im] [-2.12345+1.73512im, 0.438372-2.01187im, -0.398561+0.217356im]
Now let's prepare to solve our actual problems:
# Construct a path tracker to track the two solutions to our specific instances
# p₁ and
tracker = pathtracker(f, generic_solutions;
parameters=params,
# we just use som dummy parameters for now
targetparameters=generic_params,
startparameters=generic_params,
affine=true)
PathTracker()
Out of lack of real world data let's sample some random parameters again
specific_params = randn(length(params));
Note that the parameters order is
params
12-element Array{PolyVar{true},1}: p₁₋₁ p₂₋₁ p₃₋₁ p₁₋₂ p₂₋₂ p₃₋₂ p₁₋₃ p₂₋₃ p₃₋₃ x y z
# Setup the new parameters
set_parameters!(tracker.homotopy, generic_params, specific_params);
# track to the two solutions
r1 = track(tracker, generic_solutions[1])
r2 = track(tracker, generic_solutions[2])
PathTrackerResult{Array{Complex{Float64},1}}: • returncode → success • x → Complex{Float64}[0.454103+9.56595e-24im, 1.3872-3.60672e-23im, 0.132065+1.31272e-24im] • t → 0.0 + 0.0im • accuracy → 2.7287366859420435e-15 • accepted_steps → 17 • rejected_steps → 2
# Lets check that the tracking was successfull for both
r1.returncode == PathTrackerStatus.success && r2.returncode == PathTrackerStatus.success
true
And we get our two specific solutions:
r1.x
3-element Array{Complex{Float64},1}: -0.4432752095469262 - 1.2672493700951796e-19im -0.18955251455468036 + 9.362841922611108e-19im 0.20328610096481764 - 7.293918355685165e-20im
r2.x
3-element Array{Complex{Float64},1}: 0.45410326818290414 + 9.565951274599459e-24im 1.3871995772087389 - 3.60671834508408e-23im 0.1320653938006548 + 1.3127235321293564e-24im
We also see that these are complex vectors, but we are interested in the real solutions.
if maximum(abs ∘ imag, r1.x) < 1e-8
println("First solution: ", real.(r1.x))
end
if maximum(abs ∘ imag, r2.x) < 1e-8
println("Second solution: ", real.(r2.x))
end
First solution: [-0.443275, -0.189553, 0.203286] Second solution: [0.454103, 1.3872, 0.132065]