From acdea757d812cc2e56d46de54541635ad43b5272 Mon Sep 17 00:00:00 2001 From: Ben Montpetit <52968730+ECCCBen@users.noreply.github.com> Date: Tue, 4 Nov 2025 14:12:44 -0500 Subject: [PATCH 1/5] Changed default diagonalization_method to shur_forcedtriu in DORT --- smrt/rtsolver/dort.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/smrt/rtsolver/dort.py b/smrt/rtsolver/dort.py index f3379e51..0086aef3 100644 --- a/smrt/rtsolver/dort.py +++ b/smrt/rtsolver/dort.py @@ -143,7 +143,7 @@ def __init__( error_handling="exception", process_coherent_layers=False, prune_deep_snowpack=None, - diagonalization_method="eig", + diagonalization_method="shur_forcedtriu", ): DiscreteOrdinatesMixin.init(self, n_max_stream=n_max_stream, stream_mode=stream_mode, m_max=m_max) CoherentLayerMixin.init(self, process_coherent_layers=process_coherent_layers) From 7e5886410c567e57382a5482c4173b1c11345016 Mon Sep 17 00:00:00 2001 From: Hippolyte Signargout Date: Thu, 26 Feb 2026 15:38:29 +0100 Subject: [PATCH 2/5] Zero out non-diagonal values of diagonal npol*npol (2 or 3) blocks in addition to strictly lower triangular falues when using schur_forced_triu --- smrt/rtsolver/dort.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/smrt/rtsolver/dort.py b/smrt/rtsolver/dort.py index 0086aef3..8a8ba397 100644 --- a/smrt/rtsolver/dort.py +++ b/smrt/rtsolver/dort.py @@ -875,14 +875,22 @@ def diagonalize_shur(self, m, A, force_triu=False): # in addition it is possible to remove the 2x2 or 3x3 blocks that occurs when eigenvalues are close # forcing the lower triangular part of the shur matrix to zero solves this problem but is radical # a better algorithm would first check that the 2x2 nd 3x3 blocks are nearly diagional (values are very small) - + npol = 2 if m == 0 else 3 try: T, Z = scipy.linalg.schur(A) except scipy.linalg.LinAlgError: raise SMRTError("Schur decomposition failed.\n" + self.diagonalization_error_message()) try: if force_triu: + # Zero out strictly lower triangular values T[np.tril_indices(T.shape[0], k=-1)] = 0 + # Zero out anti-diagonal values of npol*npol diagonal blocks + if npol == 2: + T[([2 * i for i in range(T.shape[0] // 2)], [2 * i + 1 for i in range(T.shape[0] // 2)])] = 0 + else: + T[([3 * i for i in range(T.shape[0] // 3)], [3 * i + 1 for i in range(T.shape[0] // 3)])] = 0 + T[([3 * i for i in range(T.shape[0] // 3)], [3 * i + 2 for i in range(T.shape[0] // 3)])] = 0 + T[([3 * i + 1 for i in range(T.shape[0] // 3)], [3 * i + 2 for i in range(T.shape[0] // 3)])] = 0 beta, E = scipy.linalg.eig(T, overwrite_a=True) except scipy.linalg.LinAlgError: From 69af8a805f254880882dc809855d71707fffa00a Mon Sep 17 00:00:00 2001 From: Hippolyte Signargout Date: Thu, 26 Feb 2026 17:31:21 +0100 Subject: [PATCH 3/5] Correct subdiagonal zeroing in shur_forcedtriu by considering only 2*2 blocks. --- smrt/rtsolver/dort.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/smrt/rtsolver/dort.py b/smrt/rtsolver/dort.py index 8a8ba397..e53d1e8c 100644 --- a/smrt/rtsolver/dort.py +++ b/smrt/rtsolver/dort.py @@ -885,12 +885,7 @@ def diagonalize_shur(self, m, A, force_triu=False): # Zero out strictly lower triangular values T[np.tril_indices(T.shape[0], k=-1)] = 0 # Zero out anti-diagonal values of npol*npol diagonal blocks - if npol == 2: - T[([2 * i for i in range(T.shape[0] // 2)], [2 * i + 1 for i in range(T.shape[0] // 2)])] = 0 - else: - T[([3 * i for i in range(T.shape[0] // 3)], [3 * i + 1 for i in range(T.shape[0] // 3)])] = 0 - T[([3 * i for i in range(T.shape[0] // 3)], [3 * i + 2 for i in range(T.shape[0] // 3)])] = 0 - T[([3 * i + 1 for i in range(T.shape[0] // 3)], [3 * i + 2 for i in range(T.shape[0] // 3)])] = 0 + T[([2 * i for i in range(T.shape[0] // 2)], [2 * i + 1 for i in range(T.shape[0] // 2)])] = 0 beta, E = scipy.linalg.eig(T, overwrite_a=True) except scipy.linalg.LinAlgError: From 846ed07680f4ec6846a2b52ad51c68080bd30555 Mon Sep 17 00:00:00 2001 From: Hippolyte Signargout Date: Mon, 2 Mar 2026 16:20:21 +0100 Subject: [PATCH 4/5] Revert commits 7e58864 and 69af8a8 --- smrt/rtsolver/dort.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/smrt/rtsolver/dort.py b/smrt/rtsolver/dort.py index 85f58d7e..06e70cbe 100644 --- a/smrt/rtsolver/dort.py +++ b/smrt/rtsolver/dort.py @@ -983,7 +983,6 @@ def diagonalize_shur(self, m, A, force_triu=False): # in addition it is possible to remove the 2x2 or 3x3 blocks that occurs when eigenvalues are close # forcing the lower triangular part of the shur matrix to zero solves this problem but is radical # a better algorithm would first check that the 2x2 nd 3x3 blocks are nearly diagional (values are very small) - npol = 2 if m == 0 else 3 try: T, Z = scipy.linalg.schur(A) except scipy.linalg.LinAlgError: @@ -992,8 +991,6 @@ def diagonalize_shur(self, m, A, force_triu=False): if force_triu: # Zero out strictly lower triangular values T[np.tril_indices(T.shape[0], k=-1)] = 0 - # Zero out anti-diagonal values of npol*npol diagonal blocks - T[([2 * i for i in range(T.shape[0] // 2)], [2 * i + 1 for i in range(T.shape[0] // 2)])] = 0 beta, E = scipy.linalg.eig(T, overwrite_a=True) except scipy.linalg.LinAlgError: From d97e33c71d16026674ba67ad131ab29a65d0842b Mon Sep 17 00:00:00 2001 From: Hippolyte Signargout Date: Mon, 2 Mar 2026 17:17:53 +0100 Subject: [PATCH 5/5] Rename all instance of shur as schur --- smrt/rtsolver/dort.py | 42 +++++++++++++++---------------- smrt/rtsolver/test_dort.py | 6 ++--- smrt/test/test_dmrtdort.py | 6 ++--- smrt/test/test_integration_iba.py | 4 +-- 4 files changed, 29 insertions(+), 29 deletions(-) diff --git a/smrt/rtsolver/dort.py b/smrt/rtsolver/dort.py index 06e70cbe..d1e3e8de 100644 --- a/smrt/rtsolver/dort.py +++ b/smrt/rtsolver/dort.py @@ -21,17 +21,17 @@ The DORT solver is very robust in passive mode but may raise exception in active mode due to a matrix diagonalisation problem. The exception provides detailed information on how to address this issue. Two new diagonalisation approaches were added in January 2024. They are activated by setting the `diagonalization_method` - optional argument (see :py:mod:`smrt.core.make_model`). The first method (``diagonalization_method='shur'``) - replaces the scipy.linalg.eig function by a shur decomposition followed by a diagonalisation of the shur matrix. - While scipy.linalg.eig performs such a shur decomposition internally in any case, it seems that explicitly calling - the shur decomposition beforehand improves the stability. Nevertheless to really solve the problem, the second - method (``diagonalization_method='shur_forcedtriu'``) consists in removing the 2x2 and 3x3 blocks from the shur - matrix, i.e. forcing the shur matrix to be upper triangular (triu in numpy jargon = zeroing the lower part of this + optional argument (see :py:mod:`smrt.core.make_model`). The first method (``diagonalization_method='schur'``) + replaces the scipy.linalg.eig function by a schur decomposition followed by a diagonalisation of the schur matrix. + While scipy.linalg.eig performs such a schur decomposition internally in any case, it seems that explicitly calling + the schur decomposition beforehand improves the stability. Nevertheless to really solve the problem, the second + method (``diagonalization_method='schur_forcedtriu'``) consists in removing the 2x2 and 3x3 blocks from the schur + matrix, i.e. forcing the schur matrix to be upper triangular (triu in numpy jargon = zeroing the lower part of this matrix). This problem is due to the structure of the matrix to be diagonalized and the formulation of the DORT method in the polarimetric configuration. The eigenvalues come by triplets and can be very close to each other for the three H, V, U Stokes components when scattering is becoming small (or equiv. the azimuth mode 'm' is large). As a consequence of the Gershgorin theorem, this results in slightly complex eigenvalues (i.e. eigenvalues with very - small imaginary part) that comes from 2x2 or 3x3 blocks in the shur decomposition. This would not be a problem if + small imaginary part) that comes from 2x2 or 3x3 blocks in the schur decomposition. This would not be a problem if the eigenvectors were correctly estimated, but this is not the case. It is indeed difficult to find the correct orientation of eigenvectors associated to very close eigenvalues. To overcome the problem, the solution is to remove the 2x2 and 3x3 blocks. In principle, it would be safer to check that these blocks are nearly diagonal but this is @@ -119,8 +119,8 @@ class DORT(RTSolverBase, CoherentLayerMixin, DiscreteOrdinatesMixin, PlanckMixin contribution of the lowest layers is neglegible. The optical depth is a good criteria to determine this limit. A value of about 6 is recommended. Use with care, especially values lower than 6. diagonalization_method: This value set the method for the diagonalization in the eigenvalue solver. The defaut - is "eig", it uses the scipy.linalg.eig function. The "shur" replaces the scipy.linalg.eig function by a shur - decomposition followed by a diagonalisation of the shur matrix. The "shur_forcedtriu" forces the shur matrix + is "eig", it uses the scipy.linalg.eig function. The "schur" replaces the scipy.linalg.eig function by a schur + decomposition followed by a diagonalisation of the schur matrix. The "schur_forcedtriu" forces the schur matrix to be upper triangular. The "half_rank_eig" is the fastest method but requires symmetry and energy conservation which may fail with some EMModels and for some parameters. The "stamnes88" is another half rank fast method. @@ -155,7 +155,7 @@ def __init__( error_handling="exception", process_coherent_layers=False, prune_deep_snowpack=None, - diagonalization_method="shur_forcedtriu", + diagonalization_method="schur_forcedtriu", diagonalization_cache=False, rayleigh_jeans_approximation=False, ): @@ -778,8 +778,8 @@ def __init__(self, ke, ks, ft_even_phase_function, mu, weight, m_max, normalizat normalization (str or bool): Phase-function normalization option (e.g. ``'auto'``, ``'forced'``, or ``False``). symmetrization (bool): If True, enforce phase-function symmetry. - method (str): Diagonalization method. Supported values: ``'eig'``, ``'shur'``, - ``'shur_forcedtriu'``, ``'half_rank_eig'``, ``'stamnes88'``. + method (str): Diagonalization method. Supported values: ``'eig'``, ``'schur'``, + ``'schur_forcedtriu'``, ``'half_rank_eig'``, ``'stamnes88'``. cache (bool): If True, cache diagonalization results to avoid recomputation. """ @@ -808,10 +808,10 @@ def __init__(self, ke, ks, ft_even_phase_function, mu, weight, m_max, normalizat match self.method: case "eig": self.diagonalize_function = self.diagonalize_eig - case "shur": - self.diagonalize_function = partial(self.diagonalize_shur, force_triu=False) - case "shur_forcedtriu": - self.diagonalize_function = partial(self.diagonalize_shur, force_triu=True) + case "schur": + self.diagonalize_function = partial(self.diagonalize_schur, force_triu=False) + case "schur_forcedtriu": + self.diagonalize_function = partial(self.diagonalize_schur, force_triu=True) case "half_rank_eig": self.compute_half_rank_phase = True self.diagonalize_function = self.diagonalize_half_rank_eig @@ -977,11 +977,11 @@ def diagonalize_eig(self, m, A): return self.validate_eigen(beta, Eu, Ed, m) - def diagonalize_shur(self, m, A, force_triu=False): - # diagonalise the matrix. Eq (13) using Shur decomposition. This avoids some instabilities with the direct eig + def diagonalize_schur(self, m, A, force_triu=False): + # diagonalise the matrix. Eq (13) using schur decomposition. This avoids some instabilities with the direct eig # function # in addition it is possible to remove the 2x2 or 3x3 blocks that occurs when eigenvalues are close - # forcing the lower triangular part of the shur matrix to zero solves this problem but is radical + # forcing the lower triangular part of the schur matrix to zero solves this problem but is radical # a better algorithm would first check that the 2x2 nd 3x3 blocks are nearly diagional (values are very small) try: T, Z = scipy.linalg.schur(A) @@ -1256,9 +1256,9 @@ def diagonalization_error_message(self): - almost diagonal matrix. Such a matrix often arises in active mode when m_max is quite high. However it can also arises in passive mode or with low m_max. To solve this issue you can try to activate the -diagonalization_method="shur" option and if it does not work the more radical diagonalization_method="shur_forcedtriu". +diagonalization_method="schur" option and if it does not work the more radical diagonalization_method="schur_forcedtriu". These options are experimental, please report your results (both success and failure). -diagonalization_method="shur_forcedtriu" should become the default if success are reported. Alternatively you could +diagonalization_method="schur_forcedtriu" should become the default if success are reported. Alternatively you could reduce the m_max option progressively but high values of m_max give more accurate results in active mode (but tends to produce almost diagonal matrix). diff --git a/smrt/rtsolver/test_dort.py b/smrt/rtsolver/test_dort.py index 85d16791..27aad6ad 100644 --- a/smrt/rtsolver/test_dort.py +++ b/smrt/rtsolver/test_dort.py @@ -127,9 +127,9 @@ def test_shallow_snowpack(): @pytest.mark.parametrize( "microstructure_model,m_max,emmodel,diagonalization_method", - [("independent_sphere", 6, Rayleigh, "shur"), ("exponential", 16, IBA, "shur_forcedtriu")], + [("independent_sphere", 6, Rayleigh, "schur"), ("exponential", 16, IBA, "schur_forcedtriu")], ) -def test_shur_based_diagonalisation(microstructure_model, m_max, emmodel, diagonalization_method): +def test_schur_based_diagonalisation(microstructure_model, m_max, emmodel, diagonalization_method): sp = make_snowpack( thickness=[1000], microstructure_model=microstructure_model, @@ -142,7 +142,7 @@ def test_shur_based_diagonalisation(microstructure_model, m_max, emmodel, diagon nstreams = 32 # this setting fails when DORT use scipy.linalg.eig - # but this works with the shur diagonalization. Let check this: + # but this works with the schur diagonalization. Let check this: m = Model( emmodel, diff --git a/smrt/test/test_dmrtdort.py b/smrt/test/test_dmrtdort.py index ba960c12..d84bf1a4 100644 --- a/smrt/test/test_dmrtdort.py +++ b/smrt/test/test_dmrtdort.py @@ -86,8 +86,8 @@ def test_less_refringent_bottom_layer(): stickiness=0.2, substrate=make_soil("transparent", 1, 270), ) - # this test fails with some version of scipy if not using the shur method - m = make_model("dmrt_qcacp_shortrange", "dort", rtsolver_options=dict(diagonalization_method="shur_forcedtriu")) + # this test fails with some version of scipy if not using the schur method + m = make_model("dmrt_qcacp_shortrange", "dort", rtsolver_options=dict(diagonalization_method="schur_forcedtriu")) scat = active(10e9, 45) warnings.simplefilter("ignore", category=SMRTWarning) res = m.run(scat, snowpack) @@ -99,7 +99,7 @@ def test_less_refringent_bottom_layer(): def test_less_refringent_bottom_layer_VH(): # Regression test 19-03-2018: value may change if other bugs found snowpack = make_snowpack([0.2, 0.3], "sticky_hard_spheres", density=[290.0, 250.0], radius=1e-4, stickiness=0.2) - m = make_model("dmrt_qcacp_shortrange", "dort", rtsolver_options=dict(diagonalization_method="shur_forcedtriu")) + m = make_model("dmrt_qcacp_shortrange", "dort", rtsolver_options=dict(diagonalization_method="schur_forcedtriu")) scat = active(10e9, 45) warnings.simplefilter("ignore", category=SMRTWarning) res = m.run(scat, snowpack) diff --git a/smrt/test/test_integration_iba.py b/smrt/test/test_integration_iba.py index eae8f9e1..c32ff104 100644 --- a/smrt/test/test_integration_iba.py +++ b/smrt/test/test_integration_iba.py @@ -31,7 +31,7 @@ def setup_snowpack_2(): @pytest.mark.parametrize( - "method,atol", [("eig", 0.0001), ("shur", 0.0001), ("half_rank_eig", 0.0001), ("stamnes88", 0.4)] + "method,atol", [("eig", 0.0001), ("schur", 0.0001), ("half_rank_eig", 0.0001), ("stamnes88", 0.4)] ) def test_iba_dort_oneconfig_passive(setup_snowpack_2, method, atol): # prepare inputs @@ -50,7 +50,7 @@ def test_iba_dort_oneconfig_passive(setup_snowpack_2, method, atol): @pytest.mark.parametrize( - "method,atol", [("eig", 0.001), ("shur", 0.0001), ("half_rank_eig", 0.001), ("stamnes88", 1.2)] + "method,atol", [("eig", 0.001), ("schur", 0.0001), ("half_rank_eig", 0.001), ("stamnes88", 1.2)] ) def test_iba_dort_oneconfig_active(setup_snowpack_2, method, atol): # prepare inputs