Inverse Newton Method for optimization: is this the correct algorithm?

Question

I am trying to implement the algorithm in this article. I have already asked a question before about it here, and I am trying to figure out what I am doing wrong. This time, it's this section of the article that I am struggling with (bold part is mine).

[...] it can be shown that the vector $c_r$ that satisfied (29) can be estimated by minimizing the function $$\chi(c_r) = \frac{1}{2} c_r'[G(c_r) + \alpha I]c_r - c_r'\tilde{m}_r$$ where $I \in \mathcal{R}^{(s_1s_2)\times(s_1s_2)}$ is an identity matrix. The optimization for $c_r$ can be performed by using the inverse Newton method since the first and second derivative of $\chi(c_r)$ are available in closed form as $$\nabla \chi(c_r) = (G(c_r) + \alpha I)c_r - \tilde{m}_r$$ $$\nabla \nabla \chi(c_r) = G(c_r) + \alpha I$$

I looked for this inverse Newton method and didn't find anything with that exact name (not really my area of expertise). However, I did find Newton's method for higher dimensions in Wikipedia. Is that the algorithm that I should use? If not, could you give me some pointer to what I could use? I implemented it in Matlab, but my final result is still not even close to what I should have. I simulated some data, transformed it, added some noise, then tried to invert it back. I should get something like what I started with, but it's totally different and nonsensical. So I am trying to rule out possible failure areas. Here is the matlab code for all the functions, and the formula for $G(c_r)$. If anyone could double check my work, it would mean a lot. I think it looks much more difficult than it really is, but I can't wrap my head around what I'm doing wrong.

---

Formula for $G(c_r)$:

$$ G(c_r) = \tilde{K}_0 \begin{bmatrix} H(\tilde{K}'_{0,1}c_r) & 0 & \cdots & 0 \\ 0 & H(\tilde{K}'_{0,2}c_r) & \cdots & 0 \\ \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & H(\tilde{K}'_{0,N_xN_y}c_r) \end{bmatrix} \tilde{K}_0' $$

Here, $\tilde{K}_{0,i}$ is the ith column of $\tilde{K}_0$ and $H$ is the Heavyside function.

---

(Inverse?) Newton method: $x_{k+1} = x_k - \gamma [f''(x_k)]^{-1} f' (x_k)$

cr = ones(s1*s2, 1);  % Initial Guess
alpha = 1;  % Initial Guess
gamma = 1; % Step size.
prev_chi_val = 1E8;  % Just to check when it's converging
epsilon = 0.001;
maxiter = 50;

for num = 1:maxiter
   G_calc = G(cr, K0til, Nx, Ny);  % Needed for all 3 functions. Defined at the end
   chi_val = chi(cr, G_calc, alpha, mtil);  % Defined at the end
   D_chi = der1_chi(cr, G_calc, alpha, mtil); % Defined at the end
   DD_chi = der2_chi(G_calc, alpha); % Defined at the end
   cr = cr - gamma .* inv(DD_chi) * D_chi;  % Performing newton step
   fprintf('Iter %d: chi= %f, delta= %f\n', num, chi_val, prev_chi_val - chi_val);
   if abs(prev_chi_val - chi_val) < epsilon
       break
   end
   prev_chi_val = chi_val;  % Update goodness parameter
end
fprintf('Done\n')

Function to calculate $G(c_r)$

function out = G(cr, K0til, Nx, Ny)
    temp = zeros(Nx * Ny, Nx * Ny);
    for ii = 1:Nx*Ny
        temp(ii, ii) = heaviside(K0til(1:end, ii)' * cr);
    end
    out = K0til * temp * K0til';
end

Function to calculate $\chi(c_r)$

function out = chi(cr, G_calc, alpha, mtil)
    out = 0.5 * cr' * (G_calc + alpha .* eye(size(G_calc))) * cr - cr' * mtil;
end

Function to calculate $\nabla \chi(c_r)$

function out = der1_chi(cr, G_calc, alpha, mtil)
    out = (G_calc + alpha .* eye(size(G_calc))) * cr - mtil;
end

Function to calculate $\nabla \nabla\chi(c_r)$

function out = der2_chi(G_calc, alpha)
    out = G_calc + alpha .* eye(size(G_calc));
end

---

Sizes of the matrices. These are all either 2D matrices ($s_1s_2$ means multiplying these two numbers; they're the number of singular values in a previous step) or vectors, and one, $\chi$, is a number.

• $\tilde{m}_r$: ($s_1s_2 \times 1$)
• $c_r$ : ($s_1s_2 \times 1$)
• $G$ : ($s_1s_2 \times s_1s_2$)
• $\chi$ : ($1\times1$)
• $\nabla \chi$ : ($s_1s_2 \times 1$)
• $\nabla \nabla \chi$ : ($s_1s_2 \times s_1s_2$)
• $\tilde{K}_0$ : ($s_1s_2 \times N_xN_y$)|