Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

BUG: scipy 'minimize' with method='trust-constr' with equality constraints raises ValueError('expected square matrix') #20618

Closed
masinag opened this issue Apr 30, 2024 · 8 comments · Fixed by #20665
Labels
defect A clear bug or issue that prevents SciPy from being installed or used as expected scipy.optimize
Milestone

Comments

@masinag
Copy link

masinag commented Apr 30, 2024

Describe your issue.

I am trying to solve constrained optimization problems using scipy.optimize.minimize with method='trust-constr'. My problems typically involve both equality and inequality nonlinear constraints. For instance, let's say I want to solve the following problem
$\text{minimize} \quad x^4 - \frac{7}{4}x^3 - \frac{1}{2}x^2 + \frac{27}{16}x + \frac{23}{16}$
$\text{subject to} -x \leq -\frac{1}{2}, \sin(x) = 0, \cos(2x + \frac{\pi}{2}) = 0$.

I am getting errors when scipy.optimize.minimize is called, and in particolar when the QR decomposition of the Jacobian is computed within equality_constrained_sqp, since the R matrix is not squared for some reason.

This error is not raised when I only use inequality constraints (since the equality_constrained_sqp is not invoked), or when I only use one equality constraint.

I could not find much documentation apart from this link and this link. I wonder whether I am calling the function with bad parameters, in which case I think a better error should be raised, or whether this is a bug.

Reproducing Code Example

import numpy as np
import sympy as sym
from scipy.optimize import minimize, NonlinearConstraint


def jac_lambda(f, variables):
    jac_mat = sym.Matrix([[f.diff(s) for s in variables]])
    jac_lambdify = sym.lambdify([variables], jac_mat, 'numpy')
    return lambda X, *args: jac_lambdify(X)[0]


def hess_lambda(f, variables):
    hess_mat = sym.Matrix([[f.diff(s1).diff(s2) for s2 in variables] for s1 in variables])
    hess_lambdify = sym.lambdify([variables], hess_mat, 'numpy')
    return lambda X, *args: hess_lambdify(X)


def wrap_fn(expr, symbols):
    fun = sym.lambdify([symbols], expr, 'numpy')
    jac = jac_lambda(expr, symbols)
    hes = hess_lambda(expr, symbols)
    return fun, jac, hes


def main():
    f_sym = sym.sympify("x**4 + -7/4*x**3 + -1/2*x**2 + 27/16*x + 23/16")
    cc_sym = [
        sym.sympify("-x <= -1/2"),
        sym.Eq(sym.sympify("sin(x)"), 0),
        sym.Eq(sym.sympify("cos(2*x + 1/2*pi)"), 0)
    ]
    symbols = ["x"]
    constraints = []
    for c in cc_sym:
        left, right = c.args
        fun, jac, hess = wrap_fn(left, symbols)

        if c.func == sym.Eq:
            domain = (right, right)
        elif c.func == sym.LessThan:
            domain = (-sym.oo, right)
        else:
            raise ValueError("Only equality and less than constraints are supported")
        constraint = NonlinearConstraint(fun, *domain, jac=jac, hess=hess)
        constraints.append(constraint)

    fun, jac, hess = wrap_fn(f_sym, symbols)
    x0 = np.zeros(len(symbols))
    res = minimize(fun, x0, method='trust-constr', jac=jac, hess=hess, constraints=constraints)
    print(res)


if __name__ == "__main__":
    main()

Error message

Traceback (most recent call last):
  File "/mypath/examples/minimize_trust_constr.py", line 54, in <module>
    main()
  File "/mypath/examples/minimize_trust_constr.py", line 49, in main
    res = minimize(fun, x0, method='trust-constr', jac=jac, hess=hess, constraints=constraints)
  File "/mypath/venv/lib/python3.10/site-packages/scipy/optimize/_minimize.py", line 725, in minimize
    res = _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
  File "/mypath/venv/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py", line 528, in _minimize_trustregion_constr
    _, result = tr_interior_point(
  File "/mypath/venv/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/tr_interior_point.py", line 321, in tr_interior_point
    z, state = equality_constrained_sqp(
  File "/mypath/venv/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/equality_constrained_sqp.py", line 81, in equality_constrained_sqp
    Z, LS, Y = projections(A, factorization_method)
  File "/mypath/venv/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/projections.py", line 405, in projections
    Z = LinearOperator((n, n), null_space)
  File "/mypath/venv/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py", line 584, in __init__
    self._init_dtype()
  File "/mypath/venv/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py", line 182, in _init_dtype
    self.dtype = np.asarray(self.matvec(v)).dtype
  File "/mypath/venv/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py", line 236, in matvec
    y = self._matvec(x)
  File "/mypath/venv/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py", line 593, in _matvec
    return self.__matvec_impl(x)
  File "/mypath/venv/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/projections.py", line 196, in null_space
    aux2 = scipy.linalg.solve_triangular(R, aux1, lower=False)
  File "/mypath/venv/lib/python3.10/site-packages/scipy/linalg/_basic.py", line 336, in solve_triangular
    raise ValueError('expected square matrix')
ValueError: expected square matrix

SciPy/NumPy/Python version and system information

1.13.0 1.26.4 sys.version_info(major=3, minor=10, micro=12, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /usr/local/include
    lib directory: /usr/local/lib
    name: openblas
    openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=64
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.26.dev
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /usr/local/include
    lib directory: /usr/local/lib
    name: openblas
    openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= ZEN MAX_THREADS=64
    pc file directory: /usr/local/lib/pkgconfig
    version: 0.3.26.dev
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.12.0
Compilers:
  c:
    commands: cc
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  c++:
    commands: c++
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.10
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 10.2.1
  pythran:
    include directory: ../../tmp/pip-build-env-0blqy1or/overlay/lib/python3.10/site-packages/pythran
    version: 0.15.0
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: linux
Python Information:
  path: /opt/python/cp310-cp310/bin/python
  version: '3.10'
@masinag masinag added the defect A clear bug or issue that prevents SciPy from being installed or used as expected label Apr 30, 2024
@andyfaff
Copy link
Contributor

andyfaff commented May 1, 2024

The trigonometric constraints (as specified in the latex equation) are impossible to satisfy.

import matplotlib.pyplot as plt
import numpy as np

def f(x):
    return np.sin(x)

def f2(x):
    return np.sin(2*x + np.pi/2)

x = np.linspace(0, 10, 10000)

plt.plot(x, f(x))
plt.plot(x, f2(x))

plot

To be feasible the roots of f and f2 have to coincide at least one place. They will never do that because the wavelength of f is twice that of f2.

@masinag
Copy link
Author

masinag commented May 1, 2024

Actually, in my example the second constraint is cos(...) instead of sin(...), and the problem has solutions.
Figure_1

See also this discussion on StackOverflow and this smaller example by @nickodell raising the same exception. It seems that the problem is due to the presence of two equality constraints on one variable.

@masinag
Copy link
Author

masinag commented May 6, 2024

@andyfaff do you think this is indeed a bug? Do you think there is an easy way to fix it?

@andyfaff
Copy link
Contributor

andyfaff commented May 6, 2024

I'm not familiar with trust-constr, so I don't know. It'll need someone who knows.

@andyfaff
Copy link
Contributor

andyfaff commented May 6, 2024

In the meantime you could try setting breakpoints/use print statements to see how the linear operator is setup, or what the R array is on the last call

@andyfaff
Copy link
Contributor

andyfaff commented May 6, 2024

This is a peculiar edge case. If x0 has length > 1 then the issue doesn't seem to arise.

@antonior92 could you help out here?

@andyfaff
Copy link
Contributor

andyfaff commented May 7, 2024

<bear in mind that I'm not an expert in Linear Algebra>

Ok, a little bit of progress. This error arises because the number of equality constraints is greater than the number of independent variables. The problem doesn't appear if the number of equality constraints is less than or equal to the number of independent variables. I checked this by modifying a test example based on your first post - if one has two independent variables and two equality constraints there are no issues, but if there are three equality constraints there is a problem.

I then tried the example using SLSQP (1 independent variable, two equality constraints) and the minimize fails with:

message: More equality constraints than independent variables
success: False

Thus, it would not be surprising that the same behaviour is required for the trust-constr solver. I therefore don't think there's a bug.

The example can 'work' by using a different factorization_method:

minimize(f, np.ones(1), constraints=[nlc, nlc2], method='trust-constr',
         options={'factorization_method':'SVDFactorization'})

It's worth looking at the docstring there for guidance.

With respect to improving the error message that will probably require a greater knowledge of the codebase. It should be possible to catch this specific error when the projections function is called in equality_constrained_sqp. However, I'm not sure if there are other circumstances under which exceptions are raised, e.g. constraints exceeding number of independent variables for NonlinearConstraint or LinearConstraint.

@andyfaff
Copy link
Contributor

andyfaff commented May 8, 2024

How about the following. It's a little bit messy, but one can't use Exception.add_note until Python 3.11. Unfortunately it's a little fragile and depends on the error message from linalg not changing.

    try:
        Z, LS, Y = projections(A, factorization_method)
    except ValueError as e:
        if str(e) == "expected square matrix":
            # can be the case if there are more equality
            # constraints than independent variables
            raise ValueError(
                "The 'expected square matrix' error can occur if there are"
                " more equality constraints than independent variables."
                " Consider how your constraints are setup, or use"
                " factorization_method='SVDFactorization'."
            ) from e
        else:
            raise e

This would give the following error trace:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File [~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/equality_constrained_sqp.py:81](http://localhost:8889/lab/workspaces/auto-2/tree/~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/equality_constrained_sqp.py#line=80), in equality_constrained_sqp(fun_and_constr, grad_and_jac, lagr_hess, x0, fun0, grad0, constr0, jac0, stop_criteria, state, initial_penalty, initial_trust_radius, factorization_method, trust_lb, trust_ub, scaling)
     80 try:
---> 81     Z, LS, Y = projections(A, factorization_method)
     82 except ValueError as e:

File [~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/projections.py:403](http://localhost:8889/lab/workspaces/auto-2/tree/~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/projections.py#line=402), in projections(A, method, orth_tol, max_refin, tol)
    400     null_space, least_squares, row_space \
    401         = svd_factorization_projections(A, m, n, orth_tol, max_refin, tol)
--> 403 Z = LinearOperator((n, n), null_space)
    404 LS = LinearOperator((m, n), least_squares)

File [~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py:584](http://localhost:8889/lab/workspaces/auto-2/tree/~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py#line=583), in _CustomLinearOperator.__init__(self, shape, matvec, rmatvec, matmat, dtype, rmatmat)
    582 self.__matmat_impl = matmat
--> 584 self._init_dtype()

File [~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py:182](http://localhost:8889/lab/workspaces/auto-2/tree/~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py#line=181), in LinearOperator._init_dtype(self)
    181 v = np.zeros(self.shape[-1])
--> 182 self.dtype = np.asarray(self.matvec(v)).dtype

File [~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py:236](http://localhost:8889/lab/workspaces/auto-2/tree/~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py#line=235), in LinearOperator.matvec(self, x)
    234     raise ValueError('dimension mismatch')
--> 236 y = self._matvec(x)
    238 if isinstance(x, np.matrix):

File [~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py:593](http://localhost:8889/lab/workspaces/auto-2/tree/~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/sparse/linalg/_interface.py#line=592), in _CustomLinearOperator._matvec(self, x)
    592 def _matvec(self, x):
--> 593     return self.__matvec_impl(x)

File [~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/projections.py:194](http://localhost:8889/lab/workspaces/auto-2/tree/~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/projections.py#line=193), in qr_factorization_projections.<locals>.null_space(x)
    193 aux1 = Q.T.dot(x)
--> 194 aux2 = scipy.linalg.solve_triangular(R, aux1, lower=False)
    195 v = np.zeros(m)

File [~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/linalg/_basic.py:337](http://localhost:8889/lab/workspaces/auto-2/tree/~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/linalg/_basic.py#line=336), in solve_triangular(a, b, trans, lower, unit_diagonal, overwrite_b, check_finite)
    336 if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
--> 337     raise ValueError('expected square matrix')
    339 if a1.shape[0] != b1.shape[0]:

ValueError: expected square matrix

The above exception was the direct cause of the following exception:

ValueError                                Traceback (most recent call last)
Cell In[1], line 19
     13 constraints = [
     14     NonlinearConstraint(constraint1, [0, 0], [0, 0]),
     15 ]
     18 x0 = np.zeros(1)
---> 19 res = minimize(fun, x0, method='trust-constr', constraints=constraints)

File [~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/optimize/_minimize.py:725](http://localhost:8889/lab/workspaces/auto-2/tree/~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/optimize/_minimize.py#line=724), in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    722     res = _minimize_slsqp(fun, x0, args, jac, bounds,
    723                           constraints, callback=callback, **options)
    724 elif meth == 'trust-constr':
--> 725     res = _minimize_trustregion_constr(fun, x0, args, jac, hess, hessp,
    726                                        bounds, constraints,
    727                                        callback=callback, **options)
    728 elif meth == 'dogleg':
    729     res = _minimize_dogleg(fun, x0, args, jac, hess,
    730                            callback=callback, **options)

File [~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py:519](http://localhost:8889/lab/workspaces/auto-2/tree/~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/minimize_trustregion_constr.py#line=518), in _minimize_trustregion_constr(fun, x0, args, grad, hess, hessp, bounds, constraints, xtol, gtol, barrier_tol, sparse_jacobian, callback, maxiter, verbose, finite_diff_rel_step, initial_constr_penalty, initial_tr_radius, initial_barrier_parameter, initial_barrier_tolerance, factorization_method, disp)
    516         J_eq, _ = canonical.jac(x)
    517         return g, J_eq
--> 519     _, result = equality_constrained_sqp(
    520         fun_and_constr, grad_and_jac, lagrangian_hess,
    521         x0, objective.f, objective.g,
    522         c_eq0, J_eq0,
    523         stop_criteria, state,
    524         initial_constr_penalty, initial_tr_radius,
    525         factorization_method)
    527 elif method == 'tr_interior_point':
    528     _, result = tr_interior_point(
    529         objective.fun, objective.grad, lagrangian_hess,
    530         n_vars, canonical.n_ineq, canonical.n_eq,
   (...)
    538         initial_constr_penalty, initial_tr_radius,
    539         factorization_method)

File [~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/equality_constrained_sqp.py:86](http://localhost:8889/lab/workspaces/auto-2/tree/~/Documents/Andy/programming/scipy/build-install/lib/python3.10/site-packages/scipy/optimize/_trustregion_constr/equality_constrained_sqp.py#line=85), in equality_constrained_sqp(fun_and_constr, grad_and_jac, lagr_hess, x0, fun0, grad0, constr0, jac0, stop_criteria, state, initial_penalty, initial_trust_radius, factorization_method, trust_lb, trust_ub, scaling)
     82 except ValueError as e:
     83     if str(e) == "expected square matrix":
     84         # can be the case if there are more equality
     85         # constraints than independent variables
---> 86         raise ValueError(
     87             "The 'expected square matrix' error can occur if there are"
     88             " more equality constraints than independent variables."
     89             " Consider how your constraints are setup, or use"
     90             " factorization_method='SVDFactorization'."
     91         ) from e
     92     else:
     93         raise e

ValueError: The 'expected square matrix' error can occur if there are more equality constraints than independent variables. Consider how your constraints are setup, or use factorization_method='SVDFactorization'.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
defect A clear bug or issue that prevents SciPy from being installed or used as expected scipy.optimize
Projects
None yet
3 participants