Issue with CWLS in dki

[samcoveney]

I need some help figuring out if there is a problem with constrained WLS in the DKI module. See dipy/reconst/dki.py

dipy/dipy/reconst/dki.py

Lines 2611 to 2619 in e5cc120

	# Define sqrt weights as diag(yn)
	if weights:
	result = np.dot(inverse_design_matrix, y)
	W = np.diag(np.exp(np.dot(A, result)))
	A = np.dot(W, A)
	y = np.dot(W, y)
	# Solve sdp
	result = sdp.solve(A, y, check=True, solver=cvxpy_solver)

So here we see that, in order to solve a weighted least squares problem, we define W as "signal" (such that W is that square root of the weight matrix that would be applied to the squared residuals, which is "signal squared"). Then we redefine y and A so that we can solve the weight least squares problem. This is fine as a way to solve WLS, see https://en.wikipedia.org/wiki/Weighted_least_squares

The problem is that we are performing constrained weighted least squares. The positivity constraints are set up here:

dipy/dipy/reconst/dki.py

Lines 1795 to 1814 in e5cc120

	self.convexity_constraint = fit_method in {'CLS', 'CWLS', 'NLS'}
	if self.convexity_constraint:
	self.cvxpy_solver = self.kwargs.pop('cvxpy_solver', None)
	self.convexity_level = self.kwargs.pop('convexity_level', 'full')
	msg = "convexity_level must be a positive, even number, or 'full'."
	if isinstance(self.convexity_level, str):
	if self.convexity_level == 'full':
	self.sdp_constraints = load_sdp_constraints('dki')
	else:
	raise ValueError(msg)
	elif self.convexity_level < 0 or self.convexity_level % 2:
	raise ValueError(msg)
	else:
	if self.convexity_level > 4:
	msg = "Maximum convexity_level supported is 4."
	warnings.warn(msg)
	self.convexity_level = 4
	self.sdp_constraints = load_sdp_constraints(
	'dki', self.convexity_level)
	self.sdp = PositiveDefiniteLeastSquares(22, A=self.sdp_constraints)

But my understanding is that these constraints are set up with the understanding that y is log(signal). But we redefine y in the constrained weighted least squares.

As detailed in #3170 I am having problems writing the tests for CWLS (several of which were missing...) using noisy data - the tests keep failing. No problem with non-noisy data, probably because the solver detects that the problem is already solved (weighting zero residuals results in zero-residuals i.e. no effect).

For example:

    # testing return of S0 when mask is inputted
    mask_test = dki_S0F.fa > 0
    for fit_method in tested_methods:
        kwargs = {}
        if fit_method in ["CLS", "CWLS"]:
            kwargs = {"cvxpy_solver": cvxpy.CLARABEL}
        dki_S0M = dki.DiffusionKurtosisModel(gtab_2s, fit_method=fit_method,
                                             return_S0_hat=True, **kwargs)
        dki_S0F = dki_S0M.fit(DWI + np.random.normal(size=DWI.shape), mask=mask_test) # NOTE: added noise to track down the problem - if noise added, CWLS solve fails here! Weighted problem not working!
        dki_S0F_S0 = dki_S0F.model_S0

Resulting in:

_solver_opts = {}

    def solve(_solver_opts):
    
        _settings = CLARABEL.parse_solver_opts(verbose, _solver_opts)
        _solver = clarabel.DefaultSolver(P, c, A, b, cones, _settings)
>       _results = _solver.solve()
E       pyo3_runtime.PanicException: Eigval error: Eigen(1)

Note that I have rebased on master today, and the problem persists.

Given my worries about the problems with the theory (that we have supplied constraints designed for log(signal) to a solver using signal*log(signal), I think we should be worried that the constrained solver is not really working properly for this case.

edit: in addition to my theoretical worries about CWLS, note that adding noise to the signal also produces this error, regularly, for CLS... the plot thickens

[skoudoro]

Thank you for this complete report @samcoveney.

This worries me also. Can you have a look @RafaelNH and @arokem since you are more familiar with DKI?

if it is a solver issue, we might have to contact cvxpy by giving the simplest example as possible. They promoted a lot CLARABEL and deprecated ECOS.

Thank you very much.

[samcoveney]

Thanks. As a point, a lot of work has gone into #3170 which generalized the weights usable in WLS routines in both dti.py and dki.py, so if we can figure out the issue, please don't race ahead and implement a fix that it gonna be difficult for (or incompatible with) #3170 - if we can figure the problem, I can probably fold it into my PR

[samcoveney]

Sorry if not helpful, but trying to help with brain-storming the problem early on. By first doing export RUST_BACKTRACE=1 in the shell, and then running CLS on noisy signal several times until this exact fail, I can get a trace for pyo3_runtime.PanicException as the following:

thread '<unnamed>' panicked at src/solver/core/cones/psdtrianglecone.rs:464:35:
Eigval error: Eigen(1)
stack backtrace:
   0: rust_begin_unwind
             at /rustc/07dca489ac2d933c78d3c5158e3f43beefeb02ce/library/std/src/panicking.rs:645:5
   1: core::panicking::panic_fmt
             at /rustc/07dca489ac2d933c78d3c5158e3f43beefeb02ce/library/core/src/panicking.rs:72:14
   2: core::result::unwrap_failed
             at /rustc/07dca489ac2d933c78d3c5158e3f43beefeb02ce/library/core/src/result.rs:1649:5
   3: <clarabel::solver::core::cones::psdtrianglecone::PSDTriangleCone<T> as clarabel::solver::core::cones::Cone<T>>::step_length
   4: <clarabel::solver::core::cones::supportedcone::SupportedCone<T> as clarabel::solver::core::cones::Cone<T>>::step_length
   5: <clarabel::solver::core::cones::compositecone::CompositeCone<T> as clarabel::solver::core::cones::Cone<T>>::step_length
   6: <clarabel::solver::core::solver::Solver<D,V,R,K,C,I,SO,SE> as clarabel::solver::core::solver::IPSolver<T,D,V,R,K,C,I,SO,SE>>::solve
   7: clarabel::python::impl_default_py::_::<impl clarabel::python::impl_default_py::PyDefaultSolver>::__pymethod_solve__
   8: pyo3::impl_::trampoline::trampoline

and so on... (there is more to it)

Sometimes, I get failure of the solver with cvxpy.error.SolverError instead.

[samcoveney]

Okay, no movement on this. I'm chalking it up to the solver failing for now, and so i will turn off the noise in tests in #3170 that is causing failure. I will try to return to the matter in the future.

[samcoveney]

Update: I have reached out to Tom Dela Haije and Evren Özarslan, hopefully I can make some progress

Edit: @RafaelNH and @arokem any thoughts yet?