1
$\begingroup$

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

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.