1

I wrote a python code that imports some data which I then manipulate to get a nonsquare matrix A. I then used the following code to solve the matrix equation.

from scipy.optimize import lsq_linear
X = lsq_linear(A_normalized, Y, bounds=(0, np.inf), method='bvls')

As you can see I used this particular method because I require all the X coefficients to be positive. However, I realized that lsq_linear from scipy.optimize minimizes the L2 norm of AX - Y to solve the equation. I was wondering if anyone knows of an alternative to lsq_linear that solves the equation by minimizing the L1 norm instead. I have looked in the scipy documentation, but so far I haven't had any luck finding such an alternative myself.

(Note that I actually know what Y is and what I am trying to solve for is X).

Edit: After trying various suggestions from the comments section and after much frustration I finally managed to sort of make it work using cvxpy. However, there is a problem with it. First of all the elements of X are supposed to all be positive, but that is not the case. Moreover, when I multiply the matrix A_normalized with X they are not equal to Y. My code is below. Any suggestions on what I can do to fix it would be highly appreciated. (By the way my original use of lsq_linear in the code above gave me an X that satisfied A_normalized*X = Y.)

import cvxpy as cp
from cvxpy.atoms import norm
import numpy as np

Y = specificity
x = cp.Variable(22)
objective = cp.Minimize(norm((A_normalized @ x - Y), 1))
constraints = [0 <= x]
prob = cp.Problem(objective, constraints)
result = prob.solve()
print("Optimal value", result)
print("Optimal var")
X = x.value #  A numpy ndarray.
print(X)
A_normalized @ X == Y
YKZ
  • 11
  • 3
  • Have you worked out the matrix equations? You could just implement them by hand. – Mad Physicist Apr 11 '20 at 04:40
  • You might have to put that together yourself: https://stackoverflow.com/questions/49787068/l1-convex-optimization-with-equality-constraints-in-python#comment86602743_49787068 – Joe Apr 11 '20 at 05:19
  • https://stackoverflow.com/questions/58584127/l1-norm-minimization#comment103500983_58584197 – Joe Apr 11 '20 at 05:21
  • You could try this toolbox https://cvxopt.org/examples/mlbook/l1.html but I have never used it. – Joe Apr 11 '20 at 05:24
  • One solution could be to just put `return np.sum(np.dot(A_normalized, X) - Y)` in a function and wrap it in `scipy.optimize.minimize`. – Joe Apr 11 '20 at 05:28
  • From above recommendations:the (custom l1) cvxopt-solvers do not support variable-bounds as asked for and the scipy.minimize approach will lead to a non-smooth problem for some smoothness-assuming solver. Finding the right approach is highly depending on the use-case (e.g. problem.size), but cvxpy and it's available solvers (ecos for high-accuracy small-scale, scs for low-accuracy and medium-scale) are probably the easiest approach to get something which runs within a few lines of code (the code linked above can be minimized heavily as cvxpy supports `norm(x,1)`-reformulation out of the box. – sascha Apr 11 '20 at 13:26
  • After trying cvxpy it sort of works, but not quite. See the edit I made above for more details. – YKZ Apr 11 '20 at 23:40
  • Details? I don't see any and won't run your code. So i need to guess as you those details are missing. Your problem is a problem, where there exists no closed-form solution and iterative based approaches will be needed. Those iterative approaches are approaching /converging towards the optimum. They are (in theory) efficient (poly) and can be made arbitrarily accurate. In practice there is a limit and you will see, that your `>=0` is violated by something like `10^-8`. Negligible. Your mat-mul check will also have similar violations. This is normal and (in practice) will be always there. – sascha Apr 12 '20 at 09:36
  • Everyone writing `A_normalized @ X == Y` with floating-point based values, especially in (non-closed form) numerical-optimization, should really read: `What Every Computer Scientist Should Know About. Floating-Point Arithmetic`. Hint: `np.isclose()`. – sascha Apr 12 '20 at 09:40

0 Answers0