Weird behavior with PSD matrix in a Convex QP problem

回答済み

2019年12月17日 12:09 編集済み

Hi,

I am using Gurobi 9.0 to drive the computation in a manuscript about support vector machines and I found a very weird behavior. Gurobi complains that the Q matrix is not PSD (gurobipy.GurobiError: Objective Q not PSD (diagonal adjustment of 1.7e+02 would be required)).

But the spectral decomposition of the matrix using scipy.linalg says the smallest eigenvalue obtained is 5.21e-05 (positive) and the largest is 7.88 (the condition number is 151,321). So the matrix is not very ill conditioned and PSD.

Any idea of what the problem is?

I can workaround the problem by rebuilding the matrix from its decomposition (I get a matrix whose relative distance to the original one is 2.36e-14). See the commented section in the code below.

If anyone wants to download the data that exposes the error, you can grab it at

https://drive.google.com/file/d/1ooATqv6Y_UuYOEOEXaxOwbloDDYRG40y/view?usp=sharing

(I did not find a way to upload an attachment to this message)

------------ Code ------------

import gurobi as grb

import scipy as sp

import scipy.linalg as la

# Load data

data = sp.load("svm_data.npz")

C, y, K = data["C"], data["y"], data["K"]

# Verify whether the matrix is actually PSD

eigs = la.eigvalsh(K)

print("Minimum eigvalue =", eigs[0], "Maximum eigvalue =", eigs[-1], "Condition =", eigs[-1]/abs(eigs[0]))

# Create the SVM Model

n_samples = K.shape[0]

model = grb.Model("SVM")

lambd = model.addMVar(n_samples)

# Add the constraints on lambda

model.addConstr(y.reshape(1, n_samples) @ lambd ==0)

model.addConstr(lambd <= C)

# Uncomment this to workaround

# eigs, V = la.eigh(K)

# correctedK = (V * eigs) @ V.T

# print("New matrix rel error with rescpect to original =", la.norm(K - correctedK) / la.norm(K))

# K = correctedK

# Add objective

obj =0.5*(lambd @ K @ lambd) - lambd.sum()

model.setObjective(obj, grb.GRB.MINIMIZE)

model.optimize()

8件のコメント

正式なコメント
Simranjit Kaur

Gurobi Staff

2024年12月06日 00:51
This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum. Or why not try our AI Gurobot?.
Robert Luce

Gurobi Staff

2019年12月18日 21:09
Thanks for reporting, I will investigate!

Robert
0

Robert Luce

Gurobi Staff

2019年12月18日 21:18

Probably unrelated, but noteworthy nevertheless: Your matrix K contains a lot of very tiny values, e.g.,

>>> print(K)

[[ 1.00000000e+000 2.00249261e-135 -0.00000000e+000 ... 6.96007338e-032
-5.59563917e-179 3.68588049e-247]
[ 2.00249261e-135 1.00000000e+000 -1.15531076e-145 ... 3.41991649e-037
-3.76040351e-024 8.58692543e-135]
[-0.00000000e+000 -1.15531076e-145 1.00000000e+000 ... -2.10436100e-228
2.67183006e-268 -2.43865841e-029]
...
[ 6.96007338e-032 3.41991649e-037 -2.10436100e-228 ... 1.00000000e+000
-7.50509824e-074 6.20358207e-153]
[-5.59563917e-179 -3.76040351e-024 2.67183006e-268 ... -7.50509824e-074
1.00000000e+000 -9.03140556e-270]
[ 3.68588049e-247 8.58692543e-135 -2.43865841e-029 ... 6.20358207e-153
-9.03140556e-270 1.00000000e+000]]

What do these values mean? (Sometimes such artefacts point at a problem in the data generation.)

Paulo Jose da Silva e Silva

2019年12月18日 22:59
The matrix is the Gram matrix of the RBF kernel computed on the ionosphere data set after normalization of the data. Each entry represents the inner product of a nonlinear transformation of two data points. The RBF kernel as an exponential, that should help to explain the very tiny values (they should probably be zero). I believe they should be harmless.

Paulo
0
Robert Luce

Gurobi Staff

2019年12月19日 14:13 編集済み
The non-PSD error is due to a bug in gurobipy: If an arithmetic operation with an MQuadExpr object results in an FP underflow, the data underlying the MQuadExpr object may become corrupted. In your code you do
```
obj =0.5*(lambd @ K @ lambd) - lambd.sum()
```
but some of these tiny values of K times 1/2 underflow:
```
>>> np.count_nonzero(K)
970392
>>> np.count_nonzero(0.5 * K)
970140
```
A simple workaround is the clean up the tiny values of the K matrix before passing it to gurobi, or to let the underflow happen on the ndarray object before passing it to gurobipy, e.g.,
```
obj = (lambd @ (0.5*K) @ lambd) - lambd.sum() 
```
This bug will be fixed in the next maintenance release, thanks again for reporting!

Robert
0
Paulo Jose da Silva e Silva

2019年12月19日 15:06 編集済み
Hi, Bruce. Thanks for taking a look at the problem and giving some advice.

I tried to use your tips. I started trying the trick of moving the 0.5 to the obj. Actually, I wrote
```
obj = (lambd @ K @ lambd) - 2.0*lambd.sum() 
```
Which is equivalent and should also solve the problem, if I understood well. But my tests would still fail in another case. I then tried also to clean the very small values by setting to zero every entry whose absolute value is smaller than the machine epsilon times the largest entry in the matrix. It also fails in another test. In this test the matrix is only semidefinite (has 0 eigenvalues), which numerically translates in eigenvalues from -3.85907096e-14 to 6.13979525e+01. In this case Cholesky fails and, hence, maybe you are considering that this matrix should not be allowed. Anyhow, Gurobi bails out saying:

diagonal adjustment of 2.6e+03 would be required

When in fact adding the identity times 6.0e-15 already solves the problem. I uploaded the example to

https://drive.google.com/open?id=19HJppcvMckldHugFH-9SMWKRFgY5LFyL

if you want to take a look. The file name is the svm_data2.npz, you will need to edit the line that reads the file in my code above. It should fail with the 0.5 outside the quadratic and also if you move 0.5 to multiply K first.

best,

Paulo
0
Paulo Jose da Silva e Silva

2019年12月19日 15:05
Ops...

I copied the wrong link above, the correct link is

https://drive.google.com/open?id=19HJppcvMckldHugFH-9SMWKRFgY5LFyL

I will edit my post above to reflect this.

Paulo
0
Robert Luce

Gurobi Staff

2019年12月19日 20:35
Paulo,

I looked at the new data, and indeed the regularization determined by Gurobi's PSD test is far off from just a small diagonal shift. Unfortunately we run here into the limits of what our factorization based PSD test can achieve for singular matrices. As you have noticed, you are better off regularizing the matrix yourself based on spectral information before passing the optimization problem to Gurobi.

Robert
0

投稿コメントは受け付けていません。

Weird behavior with PSD matrix in a Convex QP problem

コメント

お探しのものを見つけられませんでしたか？