0

I am working to solve a given linear system of the form

$$\mathbf{A} \: \mathbf{x} = \mathbf{b}$$

where $\mathbf{x} = \begin{pmatrix} x_1 & x_2 & \cdots & x_n \end{pmatrix}^T$ with additional constraints such as

$$x_i x_j = x_k x_l \quad \text{for some } i,j,k,l \in \{1, \ldots, n \}$$ The quadratic constraints can be written in matrix form as $$\mathbf{x}^T \: \mathbf{P}_m \: \mathbf{x} = 0, \quad m \in {1, \ldots, M}$$ for say $m$ constraints. In compact form, the problem I am trying to solve is \begin{align} \widehat{\mathbf{x}} &= \arg \min \vert \vert \mathbf{A} \: \mathbf{x} - \mathbf{b} \vert \vert^2\\ &\quad \text{s.t.} \quad \mathbf{x}^T \: \mathbf{P}_m \: \mathbf{x} = 0 \end{align}

The matrix $\mathbf{P}_m$ takes the form $$\begin{pmatrix} 0 & 0 & 0 & 0 & \ldots & 0\\ \vdots & \vdots & \vdots & \vdots & \ldots & \vdots\\ 0 & \ldots & 1 & \ldots & \ldots & 0\\ 0 & 0 & 0 & 0 & \ldots & 0\\ \vdots & \vdots & \vdots & \vdots & \ldots & \vdots\\ 0 & -1 & \ldots & \ldots & \ldots & 0\\ 0 & 0 & 0 & 0 & \ldots & 0\\ \vdots & \vdots & \vdots & \vdots & \ldots & \vdots \end{pmatrix}$$ with one entry as $1$ and another as $-1$ depending upon which constraint it represents. An important point to notice, as pointed out here, is that the matrix $\mathbf{P}_m$ can be replaced by $\mathbf{P}_{m, sym} = \dfrac{\mathbf{P}_m + \mathbf{P}_m^T}{2}$ to yield the same quadratic constraint as the skew-symmetric part has no contribution. In my case, $\mathbf{P}_{m,sym}$ has both positive and negative eigenvalues.

Here are my attempts till now:

Attempt 1: Ignoring any noise in $\mathbf{b}$. I get the accurate answer with the constraints automatically satisfied without explicitly writing them, i.e. the solution is

$$\widehat{\mathbf{x}} = (\mathbf{A} \: \mathbf{A})^{\dagger} \: \mathbf{A}^T \: \mathbf{b}$$

However, as I add noise, this no longer works as the constraints are not included and the solution has more degree of freedom than is required.

Attempt 2: I tried formulating the problem as a QCQP problem which takes the following form

\begin{align} \widehat{\mathbf{x}} &= \arg \min \quad \| \mathbf{A} \: \mathbf{x} - \mathbf{b} \|^2 \\ & \quad \quad \mathbf{x}^T \: \mathbf{P}_{m,sym} \: \mathbf{x} \leq \mathbf{0} \quad \text{for } m \in \{1, \ldots, M \}\\ \end{align} for $M$ constraints. Here, I have replace the equality quadratic constraint with inequality to match the general from of the problem as given in the wiki article. As far as I understand, $\mathbf{P}_{m, sym}$ is non-convex as it has both positive and negative eigenvalues. I am trying to use this github package, built on cvxpy,to solve this non-convex problem but the solutions I get are not correct.

Attempt 3: Using method of Lagrange multipliers from Sympy package. However, it does not lead to any solution as the computation goes on forever. Here is a sample code as a minimum working example.

import sympy as sym
import numpy as np

A = np.array([[-0.968944, -0.436796, 0.196568, 0.271458, -0.430080, 0.124070], [2.150209, -0.029774, 0.019347, -0.059676, 0.008842, 0.130567], [-2.834933, 4.052180, -0.088795, -0.061756, 0.194278, -0.028226], [4.052180, 2.834933, 0.194278, -0.028226, 0.088795, 0.061756], [8.586872, -1.745449, -0.058274, 0.016471, 0.127501, 0.007528], [-1.745449, -8.586872, 0.127501, 0.007528, 0.058274, -0.016471], [-0.992538, 0.352622, 0.157250, -0.280853, -0.344055, -0.128364], [0.352622, 0.992538, -0.344055, -0.128364, -0.157250, 0.280853], [2.299443, -2.015109, 0.111921, -0.412508, -0.244876, -0.188537], [-2.015109, -2.299443, -0.244876, -0.188537, -0.111921, 0.412508], [-3.308854, 2.438692, -0.165180, 0.006687, 0.361404, 0.003056], [2.438692, 3.308854, 0.361404, 0.003056, 0.165180, -0.006687], [-0.503533, 5.245499, 0.075566, 0.620638, -0.165335, 0.283663], [5.245499, 0.503533, -0.165335, 0.283663, -0.075566, -0.620638], [4.892830, -14.009709, -0.058637, 0.312733, 0.128294, 0.142935], [-14.009709, -4.892830, 0.128294, 0.142935, 0.058637, -0.312733], [-3.866142, -3.279746, 0.221616, 0.017314, -0.484883, 0.007913], [-3.279746, 3.866142, -0.484883, 0.007913, -0.221616, -0.017314]])


b = np.array([[-0.026581], [1.695003], [0.657257], [5.048084], [4.922790], [-7.270629],[-1.589395], [1.193584], [-1.423906], [-2.574324], [-0.527349], [4.221735], [5.705906], [3.101335], [-5.224569], [-13.753907], [-5.070051], [0.165608]])


x0, x1, x2, x3, x4, x5, mu1, mu2 = sym.symbols('x0 x1 x2 x3 x4 x5 mu1 mu2', real=True)
x = sym.Matrix([[x0], [x1], [x2], [x3], [x4], [x5]])
A_sym = sym.Matrix(np.array(A))
b_sym = sym.Matrix(np.array(b))
expr = A_sym * x - b_sym
cost = sym.Transpose(expr) * expr
L = cost[0] + mu1 * (x1 * x2 - x0 * x4) + mu2 * (x1 * x3 - x0 * x5)
sol = sym.solve_poly_system([sym.diff(L, x0), sym.diff(L, x1), sym.diff(L, x2), sym.diff(L, x3), sym.diff(L, x4), sym.diff(L, x5), sym.diff(L, mu1), sym.diff(L, mu2)], x0, x1, x2, x3, x4, x5, mu1, mu2)

The solution I am expecting here is: $$x = \begin{bmatrix} 0.707107 & 0.707107 & 0.317454 & 2.633335 & 0.317454 & 2.633335 \end{bmatrix}^T$$

Question: Is the problem solvable and if yes which package can be used in Python to solve it. Any references will be appreciated!

Zero
  • 464
  • This may be helpful, if you have enough points, and the noise level is relatively low. One of the answers even includes efficient Numpy code https://math.stackexchange.com/q/99299/207316 However, it doesn't address the issue of additional constraints. – PM 2Ring Sep 06 '21 at 10:53
  • @RodrigodeAzevedo Hi, it means all eigenvalues are zero just like a zero matrix with all eigenvalues $0$ which makes it both negative semidefinite and positive semidefinite...it is the same logic as saying $0$ is both non-negative and non-positive – Zero Sep 06 '21 at 12:24
  • @Zero You should not pay attention to the eigenvalues of $P_i$ if $P_i$ is not symmetric. The skew-symmetric part contributes $0$. Focus on the symmetric part. – Rodrigo de Azevedo Sep 06 '21 at 12:27
  • @RodrigodeAzevedo the reason for mentioning the eigenvalues of $\mathbf{P}_i$ was to point out that the problem is convex as stated here: https://en.wikipedia.org/wiki/Quadratically_constrained_quadratic_program – Zero Sep 08 '21 at 11:53
  • @PM2Ring Yeah I am looking to see if there is a formulation in cvxpy that can solve this with constraints... – Zero Sep 08 '21 at 11:53
  • @Zero Again, your approach is misguided. If your eigenvalues are not of a symmetric matrix, you should dismiss them. Again, focus on the symmetric part. – Rodrigo de Azevedo Sep 08 '21 at 12:05
  • @RodrigodeAzevedo I might be out of depth but what symmetric part are you referring to? It is a linear equation with quadratic constraints...if I do not consider the constraints dictated by $\mathbf{P}_i$, all I have is a least squares solution which will give me the exact solution satisfying the constraint when no noise is present but would fail to do so in case of noise. So can you elaborate which symmetric part you are referring to because I have no idea...my apologies if this should have been straight-forward and I cannot get it somehow :) – Zero Sep 09 '21 at 12:46
  • @Zero Take a look at this. – Rodrigo de Azevedo Sep 09 '21 at 12:57
  • @Zero Since no one answered the question, there is no need to append an edit. Just rewrite the whole question and orphan this comment section. How did you go from equality constraints to inequality constraints? Have you considered using Lagrange multipliers and symbolic (rather than numerical) methods? How large is $n$? – Rodrigo de Azevedo Sep 09 '21 at 17:26
  • @RodrigodeAzevedo Thanks for bearing with me. I have edited the question to include all I know. The inequality is just to write it into a general form as I have seen in texts and the wiki page. In my case, $n = 6$, but it can go higher, which probably will make it impractical to solve symbolically. I can use Lagrange multipliers but the method I know will only give a local minima and would depend on the initialization. – Zero Sep 09 '21 at 21:24
  • @Zero It may be possible to use symbolic methods. Take a look at this. One can also use convex relaxations and try to solve it numerically. Take a look at this. – Rodrigo de Azevedo Sep 09 '21 at 21:37
  • @RodrigodeAzevedo I have tried the method of Lagrange multipliers but the computation seems to go on forever. I have added the minimum working example if you want to take a look...thanks for your comments till now :) – Zero Sep 12 '21 at 18:46
  • @Zero Interestingly, you used NumPy arrays instead of SymPy matrices. I had never encountered SymPy code with NumPy arrays. I am not claiming that had your code used SymPy matrices it would have produced the desired results before the Sun turns into a red giant. – Rodrigo de Azevedo Sep 12 '21 at 19:04
  • @RodrigodeAzevedo Well I quickly tried right now assuming that was the problem. But unfortunately that does not help :) – Zero Sep 12 '21 at 19:19
  • @Zero Do you have access to Mathematica? When working with systems of polynomial equations, I usually use Macaulay2. – Rodrigo de Azevedo Sep 12 '21 at 19:24
  • @RodrigodeAzevedo I do not have access to Mathematica (unless I include the version I have in my Raspberry but I am afraid it wouldn't be as powerful)...I do have access to Maple and thus Macaulay2 sounds interesting...If I find a solution or at least why I cannot find the correct solution, I'll update it here...thank you for all your input till now :) – Zero Sep 12 '21 at 20:35
  • @Zero If it does not work, one can try to relax. – Rodrigo de Azevedo Sep 12 '21 at 20:36

0 Answers0