Default DORT diagonalization_method for active sensors

[ECCCBen]

When simulating a wide range of snow conditions in an active setup, some snow profiles do not converge to a solution, thus generating NaNs.

This can be reproduced with the codes/data available here: TVCExp18-19.

The rtsolver_option = {'diagonalization_method'='shur_forcedtriu'} option fixed the issue without changing the results of the study TVCExp18-19/diag

This option was also used in recent MCMC SWE retrieval using airborne Ku SAR data: TVCExp18-19_SWE.

I suggest making shur_forcedtriu the default option.

[ghislainp]

Thank you for reporting success with this option.

To complement: the rational behind shur_foredtriu is explained in a note in the documentation: https://smrt.readthedocs.io/en/v1.5.1/smrt.rtsolver.html#module-smrt.rtsolver.dort .
It was implemented nearly 2 years ago.

The reason why it is not the default yet is:

• as indicated in the note, I'm waiting for several reports of success. In fact I'd like not only reports that it successfully removes NaNs, but also reports that it provides the same results as the eig default (robust) method when both are working well, on a wide range of cases.
• this options is a dirty hack, applying a brute force patch that solves the problem in some identified conditions (see the note), but I'm not sure it does not degrade the results when applied in other conditions that also produces Nan (but for other unidentified conditions) . The reason I'm not confident is because I've well identified the problem (I believe) -- this is due to the Gershgorin theorem -- but I've never found a proper solution, neither in the literature on RT solvers (Stamnes, etc...), nor in the math literature on this theorem.

My suggestion to move on is to get help from mathematicians to look at the problem (I can describe it, it is simple). They might find a better and reliable solution, that we can trust immediately. If not, as physicists, we'll rely on many simulations that give confidence (in Nans and non-Nans conditions)...

In brief, I mean that we need a bit more before setting this option as default, but I encourage the community to perform tests with both the eig and shur_foredtriu options on their use case, and compare the results.

[ECCCBen]

This is reasonable.

As I mentioned briefly, I optimized the polydispersity (Picard et al., 2022) using a least square method (Montpetit et al., 2024), and some values of polydispersity created NaNs with eig. I had developed a workaround initially to converge to proper polydispersity values, but know with the shur_forcedtriu option, I converge to the same values without the need for the NaNs workaround.

For my MCMC SWE retrieval work (Montpetit et al., 2025), it was not possible to use SMRT with eig since the method implemented with pyMC would simply fail. Setting the method to shur_forcedtriu fixed this issue and made it possible to retrieve not only SWE but realistic snow profiles, comparable to measurements in the field, with SMRT at the core of the MCMC model.

As you said, I will keep testing both in my work, and I encourage others to try it out.

[hsignarg]

Multiple tests have been carried out with shur_forcedtriu and show no sign of error.
Additionally, I checked the maximal theoretical error coming from zeroing out unintentional non-zero values. It turns out to be similar to unavoidable machine error.

The math behind it is the following: the Schur triangulation in the diagonalisation algorithm outputs 2*2 blocks corresponding to close conjugate complex eigenvalues when the real input matrix is not diagonalisable. When the matrix is diagonalisable in the real space but two eigenvalues are very close, it may happen that the algorithm outputs complex conjugate eigenvalues with imaginary part close to machine error. The error on the real part of the eigenvalues is of the same order.

While keeping these imaginary parts can yield errors (or nans), removing it by zeroing out non-diagonal entries of 2*2 blocks is only source of unsignificant error on the real eigenvalues.

This is not what is actually done in shur_forcedtriu but it is what will be done with #71

[hsignarg]

It turns out after testing the method I proposed in #71 that I slightly misunderstood the output of the Schur algorithm in the case of close eigenvalues.

As I mentioned, it may happen that when two eigenvalues $\lambda_0$ and $\lambda_1$ are close, instead of returning a $2\times 2$ diagonal block

$$\begin{bmatrix}\lambda_0&\beta\\0&\lambda_1\end{bmatrix}$$

the Schur algorithm outputs a slightly incorrect

$$\begin{bmatrix} \alpha & \beta \\ \gamma & \alpha \end{bmatrix}$$

representing two complex conjugate eigenvalues with real part $\alpha$.

However, I misunderstood what their imaginary part was: it is actually $\pm\sqrt{-\gamma\beta}$, which means that even though both $\gamma$ is necessarily insignificant and can be zeroed out and $\sqrt{-\gamma\beta}$, which is also the error on the eigenvalues, is small, $\beta$ can actually be quite large and is not $0$ in the case of real eigenvalues. It should thus not be zeroed out.

This means the initial forcedtriu method of #71 was actually correct.