5
$\begingroup$

(I have asked this question on StackOverflow previously but it has been pointed to me that CSSE or MSE could be more appropriate)

I have to solve a constrained optimization problem of the following form, where the only unknown is $x$:

$$ x^{*} = \arg \min_{x} \left \| Ax - b \right \|^{2} \qquad x \in \mathbb{R}_{\geq 0}^{n} , ~ A \in \mathbb{R}_{\geq 0}^{n \times n} , b \in \mathbb{R}_{\geq 0}^{n} $$

in other words a non-negative least squares problem (NNLS). Alternatively, I could solve a linear program (aware that these are not the same, but I would like a solution from either, whichever is more convenient):

$$ \text{minimize} \sum_i \xi_i $$ $$ \text{subject to:} ~~ Ax = b ~ + \xi$$ $$ x \in \mathbb{R}_{\geq 0}^{n} , \xi \in \mathbb{R}_{\geq 0}^{n}$$

Now so far so good. My problem is that the matrix A and the vector b that I am using contain extremely small entries (1e-60,1e-100). Note that all numbers are roughly this small. This is because they come from evaluations of a high dimensional pdf. As far as I can tell, not even the most precise solvers could properly handle such numbers. The rest of my algorithm deals with these numbers just fine, since all operations are carried out in log-space, as it is commonplace with probabilities.

Trying to solve the original problem with either method I presented, that is using e.g. scipy.optimize.nnls or scipy.optimize.linprog, results in the solver simply returning a vector of zeros.

One could think about solving the following modified problem (for example):

$$ x^{**} = \arg \min_{x} \left \| \log (A) x - \log(b) \right \|^{2} \qquad x \in \mathbb{R}_{\geq 0}^{n} , ~ A \in \mathbb{R}_{\geq 0}^{n \times n} , b \in \mathbb{R}_{\geq 0}^{n} $$

An analogue modification could be done for the LP showed earlier. While this does not run into the same optimization issues, the optimal solution to this modified problem is not the same as one for the original problem. That is, $x^{*} \neq x^{**} $ and also $x^{*} \neq \exp(x^{**}) $. Solving this modified problem and exponentiating its solution does not give completely nonsensical results, but it is not good enough for my purposes.

How would I go about solving the original problem despite the optimization issues given by the small entries of $A$ and $b$?

$\endgroup$
  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$ – nicoguaro Jun 19 at 15:05
  • 6
    $\begingroup$ There is no problem with having small numbers -- just multiply both $A$ and $b$ by $10^{16}$ and you are dealing with numbers around 1. There is only an issue if you have both small and large numbers. $\endgroup$ – Wolfgang Bangerth Jun 20 at 23:16
  • 2
    $\begingroup$ Could you make an example problem available (e.g. by providing a link to files containing A and b)? $\endgroup$ – Brian Borchers Jun 24 at 0:23
0
$\begingroup$

To tackle your issue of very, very small numbers, you need to use a arbitrary precision arithmetic library, like MPFR. https://www.mpfr.org/

MPFR is awesome and will continually jack up the precision until it is sufficient to avoid roundoff error or you run out of memory. In my experience, I never used more than a 128 bit mantissa (a double has something like 53). Goodbye numerical limitations! Your program will run slower, but it will succeed.

If your preferred solver does not support this data type, you can write your own with a very simple gradient descent implementation and a change of variables.

Let $x_i=y_i^2$. The vector function is now expressed as $F=Ay^2-b$. The gradient descent iteration goes something like this:

$$y_{new}=y_{old}-\gamma∇||F||^2$$ where $∇||F||^2=2(∇F)^TF$ and $∇F_{ij}=2A_{ij}y_j$ and $\gamma$ is an arbitrary (usually small) positive step size. Recover $x$ by squaring the values of $y$, which is always non-negative. There are sophisticated ways of varying $\gamma$ each step to improve convergence. No matter what you do, you will most likely have to try a few (or many) random starting points to find the best solution.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.