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

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$$)|