Possibly wrong derivative in `VectorAmplitude` about simpeg

[lheagy]

@johnweis0480 has run into some questions here to. Would you have an example to share?

[johnweis0480]

I was looking at this a while back and saw that the derivative test for the smoothness would fail.

What I initially implemented and suggested last spring was to take the amplitude of the cell gradients of the components of magnetization instead of taking the cell gradient of the amplitude. I believe that the failure is because the derivative is actually matching my initial implementation of this regularization for the smoothness.

Currently, the f_m of the vector smoothness class is doing this:

Instead, I believe is should be doing this

to match the derivative and my initial implementation. By changing to the second, the derivative test will pass.

I do not believe that the format of f_m is well defined as it provides no information on the direction of magnetization. I did compare using cartesian smoothness (top row) to my implementation (bottom row) for MVI with positive results.

[jcapriot]

@johnweis0480 What you're suggesting is taking the amplitude of the gradient of the vector components, which is not the same thing as the gradient of the vector amplitudes.

[domfournier]

This is not really related to what @thibaut-kobold is raising, but the result looks good. You might be onto something. I think we should write down the math for this.

@jcapriot Are we not running derivative tests for the regularization? It should brake if @johnweis0480 is right.

[jcapriot]

I thought we were, but I'll need to double check. What @johnweis0480 is proposing is a JTV measure on the vector components of the model (Which is a valid thing to do to encourage joint sparsity in the vector components, just not explicitly what I think the VectorAmplitude regularizer is trying to do.

[lheagy]

Hi all, @jcapriot, you are correct in that @johnweis0480 is suggesting taking the amplitude of the gradient of the vector components, and I support us doing this as I think it makes more sense than looking at spatial derivatives of the amplitude. I will expand…

The current implementation for the objective function is

This has some potential problems. Let’s consider 2 cells right next to each other. Cell 1: [1, 0, 0], cell 2: [0, 1, 0]. If we take the amplitude, and then spatial derivatives, the value of the objective function is 0, so we have reached a minimum. But I don’t think that is the behaviour that we want.

Instead, if we take the amplitude of the vector differences as John has suggested,

then in the example, I pose Cell 1: [1, 0, 0], cell 2: [0, 1, 0], the value of the objective function would be 2, and there would still be progress to be made in the minimization.

So the latter approach would promote smoothly varying vectors, which I think is the behaviour we want.

What makes this a bit subtle is that the derivatives as implemented are not correct, (@johnweis0480: do you have an example that you could share that runs the deriv test?). I don’t see a derivative test in the current SimPEG test suite (I checked in test_vector_amplitudes.py and in test_regularization

The derivate as implemented is first taking gradients

        return sp.block_diag([self.cell_gradient] * self.n_comp) @ self.mapping.deriv(
            self._delta_m(m)
        )

https://github.com/simpeg/simpeg/blob/211f627874c46d3a03450b840970b9a10a7fbe4b/SimPEG/regularization/vector.py#L947C37-L949

and then amplitudes

        return self.f_m_deriv(m).T * (
            self.W.T
            @ self.W
            @ (self.f_m_deriv(m) @ d_m).reshape((-1, self.n_comp), order="F")
        ).flatten(order="F")

https://github.com/simpeg/simpeg/blob/211f627874c46d3a03450b840970b9a10a7fbe4b/SimPEG/regularization/vector.py#L426C8-L430C29

So the direction that the gradients are taking us is more aligned with the approach of taking spatial gradients and then amplitudes.

[jcapriot]

I don't mean to imply that @johnweis0480's method is not a good thing, my point is it's not what I believe this class was meant to do, say if you only have reference information on the magnitude of the vector, the approach of taking the amplitude of the gradient of the vector components will not correspond to that information.

[johnweis0480]

plot_inv_mag_MVI_VectorAmplitude_failed_test.zip

For reference, here is a copy of the file used for the test and a printout of the result.

[domfournier]

simpeg/SimPEG/regularization/vector.py

Line 429 in 211f627

@ (self.f_m_deriv(m) @ d_m).reshape((-1, self.n_comp), order="F")

@domfournier : For the first derivative deriv, we are calling twice f_m_deriv. I think the second one is potentially wrong. It should be f_m (and potentially need then to remove d_m)

See for reference the derivative of the base class: https://github.com/simpeg/simpeg/blob/211f627874c46d3a03450b840970b9a10a7fbe4b/SimPEG/regularization/base.py#L502C17-L502C17

Back to @thibaut-kobold original question - After reviewing, the algebra seems to work, but the confusion comes from what's in f_m_deriv. We are currently mixing the purpose of the derivative of the mapping, and derivatives of the kernel function. We are usually using the mapping to handle wires (grabs a chunk of the model vector for the regularization), and assume that this matches the shape of the kernel function f_m, which is not the case here. @jcapriot and I managed to re-use the mechanics of the base class and reduce code duplication, but it is less clear. I will submit a PR shortly to make it more explicit.