r-barnes / faster-unmixer

A faster implementation of the sediment inverse/unmixing scheme proposed in Lipp et al (2021).
5 stars 1 forks source link

Fix DPP #21

Closed r-barnes closed 1 year ago

r-barnes commented 1 year ago

I've explained this upstream here and here. Copying from the cvxpy repo:

In the following test program

import cvxpy as cp

def cp_log_ratio_norm(a, b):
  # Both `a * cp.inv_pos(b)` and `a / b` make this problem non-DPP
  return cp.maximum(a * b, b * cp.inv_pos(a))

var = cp.Variable(pos=True)
param = cp.Parameter(pos=True)
param.value = 5
objective = cp.Minimize(cp_log_ratio_norm(var, param))
problem = cp.Problem(objective, [])
objective_value = problem.solve()

print(f"Objective value = {objective_value}")
print(f"Status = {problem.status}")

we see that using the reciprocal of a parameter causes DPP to break with the message:

miniconda3/envs/faster-unmixer/lib/python3.10/site-packages/cvxpy/reductions/solvers/solving_chain.py:178: UserWarning: You are solving a parameterized problem that is not DPP. Because the problem is not DPP, subsequent solves will not be faster than the first one. For more information, see the documentation on Discplined Parametrized Programming, at
    https://www.cvxpy.org/tutorial/advanced/index.html#disciplined-parametrized-programming

Looking at the docs we find that this is expected:

As another example, the quotient expr / p is not DPP-compliant when p is a parameter, but this can be rewritten as expr * p_tilde, where p_tilde is a parameter that represents 1/p.

But what if we need both p and 1 / p? Keeping those two synchronized is challenging. To fix this, we can build a class ReciprocalParameter that introduces parameters for both p and 1 / p and keeps them in sync, like so:

import cvxpy as cp
from typing import Optional

class ReciprocalParameter:
    """Used for times when you want a cvxpy Parameter and its ratio"""

    def __init__(self, *args, **kwargs) -> None:
        self._p = cp.Parameter(*args, **kwargs)
        # Reciprocal of the above
        self._rp = cp.Parameter(*args, **kwargs)

    @property
    def value(self) -> Optional[float]:
        """Return the value of the Parameter"""
        return self._p.value

    @value.setter
    def value(self, val: Optional[float]) -> None:
        """
        Simultaneously set the value of the Parameter (given by `p`)
        and its reciprocal (given by `rp`)
        """
        self._p.value = val
        self._rp.value = 1 / val if val is not None else None

    @property
    def p(self) -> cp.Parameter:
        """Returns the parameter"""
        return self._p

    @property
    def rp(self) -> cp.Parameter:
        """Returns the reciprocal of the parameter"""
        return self._rp

def cp_log_ratio_norm(a, b: ReciprocalParameter):
  # Both `a * cp.inv_pos(b)` and `a / b` make this problem non-DPP
  return cp.maximum(a * b.rp, b.p * cp.inv_pos(a))

var = cp.Variable(pos=True)
param = ReciprocalParameter(pos=True)
param.value = 5
objective = cp.Minimize(cp_log_ratio_norm(var, param))
problem = cp.Problem(objective, [])
objective_value = problem.solve()

print(f"Objective value = {objective_value}")
print(f"Status = {problem.status}")