How to construct constraints code?

[ericzzj1989]

Great work!

I am looking at the code here

usbs/utils/qap_helpers.py

Line 195 in f36da24

constraint_indices = b.shape[0] + jnp.tile(jnp.arange(0, n-1)[:, None], (1, 3)).reshape(-1,)

I am wondering if you could provide some clarification on how the constraints in the QAP are constructed. Specifically, could you explain the logic or steps involved in defining the constraints at that point in the code?

Thanks!

[rangell]

That block of code is part of the vectorized implementation of the constraint "diag(Y) = vec(B)" from the paper. The more intuitive code (in the for loop) is in the comment above that block of code. Basically, this code is just transforming the constraints -- in this case, diag(Y) = vec(B) -- into the matrices for the linear operator that USBS requires to run.

[ericzzj1989]

Thanks for your previous help!

I’ve successfully implemented my sdp_state, initialize_state function, and custom helpers.py script, and everything runs without errors.

However, I’m wondering how you typically validate the correctness of the SDP optimization. Specifically, what specific printed information here should I look for to ensure the optimization is progressing correctly?

Thanks in advance for your insights!

[rangell]

Great!

Generally, obj_gap, infers_gap, and max_infeas are used to determine the optimization is progressing correctly. An optimal solution has been found when all of these are 0. Often times we just stop optimization when these are "close enough" to 0.

Let me know if you have any further questions.

[ericzzj1989]

Thanks for the helpful explanation!

I also noticed that here, the optimization progress is monitored using obj_gap, infers_gap, and max_infeas. I ran my implementation and observed that these three values decrease as the program runs, indicating that the optimization seems progressing as expected.

However, I’ve encountered a question regarding the final output. I've fixed the input for my optimization problem, but when the program finishes running, the resulting sdp_state.P seems to vary each time I run it, even though the optimization seems to converge (i.e., obj_gap, infers_gap, and max_infeas are all decreasing).

I wonder if it’s correct to view sdp_state.P as the final optimized parameter, or if there might be some other considerations I should take into account when interpreting the results?

Thanks again for your insights!

[rangell]

The sdp_state.P parameter is a projection of the primal variable (i.e. "X") using the sketching matrix sdp_state.Omega. You can use the reconstruct function in utils/common.py to get a low-rank factorized approximation of the primal solution.

There are several reasons why sdp_state.P might be different every run. One explanation is sdp_state.Omega is different -- sdp_state.Omega is randomly generated, and hence depends on the random seed. Another possibility is that there are multiple solutions to your SDP. If the difference is fairly small, this could be due to even numerical error that we don't have control over.

[ericzzj1989]

Thanks again for your useful explanation regarding the randomness, especially the part about sdp_state.Omega being randomly generated. That clarifies a lot about the variability in the results.

However, I’ve encountered another issue when printing sdp_state.X. It shows None — I’m wondering why this might be happening?

To provide more context, my optimization problem has a C matrix of size 16x16, so the optimization parameter is 16-dimensional. After the program finishes running, I printed the X, P, and y values from the sdp_state, and here are the results:

Final X: None

Final P:
[[-3.56997478e+00  4.47096403e+00 -1.91077790e+00  6.72280013e+00
  -3.05876587e+00 -2.08602868e-01  5.89763859e-01  1.66677713e+00
   7.45438122e+00 -9.70609889e+00  1.80906283e+00  3.88297424e+00
   4.94837924e+00  2.93500055e+00 -4.36636081e+00 -2.74440901e+00]
 [-4.78898292e+00  6.00717057e+00 -2.57621014e+00  9.03676007e+00
  -4.11318694e+00 -2.71287480e-01  7.91740612e-01  2.22722684e+00
   1.00282743e+01 -1.30534205e+01  2.43961855e+00  5.23109136e+00
   6.64274589e+00  3.93443930e+00 -5.86623582e+00 -3.68487304e+00]
 [ 2.66585937e+00 -3.34615561e+00  1.43656194e+00 -5.03237175e+00
   2.29179927e+00  1.54165653e-01 -4.37577836e-01 -1.24058820e+00
  -5.58401740e+00  7.26555761e+00 -1.35989041e+00 -2.91302464e+00
  -3.69888925e+00 -2.19364948e+00  3.27073242e+00  2.05530316e+00]
 [-4.49960373e+00  1.07738822e+01 -8.60828060e+00  1.55433946e+01
  -9.19518536e+00 -1.54096834e+00 -2.86627385e+00  1.12197314e+00
   1.90867631e+01 -2.00648981e+01  7.65166559e+00  1.18934910e+01
   8.81959590e+00  7.12722678e+00 -1.38491196e+01 -9.04695243e+00]
 [-9.66795077e-01  1.05161725e+01 -1.17400256e+01  1.46142446e+01
  -1.04942010e+01 -2.38106692e+00 -6.38405676e+00 -1.31683352e+00
   1.95539094e+01 -1.67894164e+01  1.02122568e+01  1.37180017e+01
   6.01693436e+00  7.00896285e+00 -1.62980548e+01 -1.08650972e+01]
 [ 1.30938293e+00 -4.22667029e+00  3.82092938e+00 -6.02205649e+00
   3.80977580e+00  7.24677738e-01  1.60411094e+00 -1.21703823e-01
  -7.60780551e+00  7.49376204e+00 -3.36597377e+00 -4.94352487e+00
  -3.11509685e+00 -2.80603276e+00  5.80528851e+00  3.82235869e+00]
 [ 4.31297718e+00  2.69248172e+00 -7.45257416e+00  3.01579199e+00
  -4.70602521e+00 -1.73786494e+00 -6.38319529e+00 -3.57625178e+00
   6.27034222e+00 -6.40898312e-01  6.27729998e+00  6.35423529e+00
  -1.91034522e+00  1.82988059e+00 -7.84858272e+00 -5.46250827e+00]
 [ 1.73067711e+00  1.08169822e+00 -2.99217588e+00  1.21239257e+00
  -1.88890505e+00 -6.96576969e-01 -2.56125671e+00 -1.43739795e+00
   2.51898617e+00 -2.61161206e-01  2.52098201e+00  2.55313582e+00
  -7.65873319e-01  7.33769323e-01 -3.15071227e+00 -2.19220136e+00]
 [-1.38418803e+00 -8.63133271e-01  2.39071915e+00 -9.67053790e-01
   1.50902562e+00  5.56193936e-01  2.04696901e+00  1.14876331e+00
  -2.01153034e+00  2.06415825e-01 -2.01398209e+00 -2.03924245e+00
   6.13608143e-01 -5.85344585e-01  2.51639567e+00  1.75102458e+00]
 [-2.81139590e+00  9.30392760e+00 -8.48135611e+00  1.32495797e+01
  -8.41841506e+00 -1.60260348e+00 -3.59847451e+00  2.12355092e-01
   1.67765842e+01 -1.64558315e+01  7.46836982e+00  1.09341433e+01
   6.80316718e+00  6.16911331e+00 -1.28298983e+01 -8.44947276e+00]
 [-5.79538543e+00  6.55265764e+00 -2.25250783e+00  9.94616086e+00
  -4.23534136e+00 -1.60665539e-01  1.45429573e+00  2.84130591e+00
   1.07760213e+01 -1.46868217e+01  2.20037741e+00  5.34385165e+00
   7.67850573e+00  4.29345126e+00 -5.94273931e+00 -3.68857639e+00]
 [ 1.29285389e+00 -4.22088285e+00  3.83099074e+00 -6.01506318e+00
   3.81102453e+00  7.19727401e-01  1.61225358e+00 -1.06050961e-01
  -7.60865003e+00  7.48442629e+00 -3.37519409e+00 -4.95248706e+00
  -3.09961260e+00 -2.79591294e+00  5.80454947e+00  3.82076494e+00]
 [-8.36306055e-02  1.56733700e-02  5.99123789e-02  4.53784316e-02
   2.32239084e-02  4.12785826e-02  9.27502825e-02  3.76378596e-02
   2.79112585e-02 -1.32497736e-01 -4.43256002e-02 -8.37788847e-03
   7.06618863e-02 -1.25516794e-02  5.80663423e-02  4.95709959e-02]
 [-2.04348762e-02  3.25097901e-03  1.56449718e-02  1.11125380e-02
   7.53348443e-03  1.40003610e-02  2.56341331e-02  5.89656098e-03
   8.00125861e-03 -3.72087195e-02 -1.07308403e-02  6.26780704e-04
   1.69666184e-02 -7.20234478e-03  1.82453538e-02  1.61824620e-02]
 [-1.60115996e-02  1.53160260e-03  1.45108125e-02  3.19894504e-04
   4.50938355e-03 -6.94827279e-03  8.49817951e-03  1.87113868e-02
  -1.05801252e-02  2.06654118e-03 -1.39078468e-02 -1.70911753e-02
   1.29033168e-02  1.01104208e-02  3.45118035e-03  1.31361235e-04]
 [ 1.29045663e+00 -4.22092236e+00  3.83287024e+00 -6.01440805e+00
   3.81134819e+00  7.19395108e-01  1.61413940e+00 -1.03432399e-01
  -7.60937449e+00  7.48361797e+00 -3.37704783e+00 -4.95476817e+00
  -3.09781483e+00 -2.79484062e+00  5.80537101e+00  3.82096468e+00]]

Final y:
[-0.05819907  0.01262477  0.10890626 -0.0764923   0.21819564 -0.01922856  0.07592298  0.20026329  0.04454488  0.05958619 -0.39044204 -0.19407816  0.94357158 -2.04478542 -0.74668516  3.39870364 -0.160199    0.32517152  3.37454917 -1.35958547 -0.54030346  1.57413383 -1.174989   -0.36173332  1.63159657  0.          7.74398061]

From the output, I noticed that the Final P matrix has a dimension of 16x16x16. I’m curious about why the output is structured like this and how I should interpret it. Is this the expected result for P?

Also, do you have any suggestions on how I can further inspect the results, especially for X, or if there’s something I might be missing in my implementation?

Thanks for your time and support!

[rangell]

If your problem size allows you to store the entire primal matrix sdp_state.X, you can turn off sketching by setting --sketch_dim=0. Turning off sketching will make sdp_state.P equal to None.

If your problem is bigger, then you have to turn sketching back on (setting sketch_dim to some number greater than 0, higher being more accurate). In order to inspect the final answer, you need to reconstruct_from_sketch to get a low-rank approximation of the optimal solution.

I think for QAP, we know that the solution has a specific rank, so we just set that permanently and round using Munkres. But, you can easily turn this off if you want to inspect the entire primal matrix.

Let me know if you have any further questions.

[ericzzj1989]

Thanks again for your previous response! It was really helpful.

Right now, I’m working on accelerating the USBS optimization speed, specifically aiming to speed up each iteration. I noticed that reducing the rho value helps improve the iteration speed, and I’m also experimenting with adjusting other parameters.

Do you have any recommendations on which parameters would be most effective for further improving iteration speed? Or any general suggestions on optimizing the performance of each step?

Looking forward to your insights!

[ericzzj1989]

Thanks again for your previous response!

I want to follow up with a new question based on my latest progress.

Initially, I wasn’t using sketching, which led to slower optimization. Now that I’ve incorporated sketching, I’ve seen a speed improvement. However, this has introduced a new question regarding the final output.

Since I’m using sdp_state.P as the final result instead of sdp_state.X, I noticed that performing an SVD decomposition on sdp_state.P only produces sketch_dim singular values, rather than the original n dimensions.

How can I project resulting sdp_state.P back to the original n-dimensional parameter space? Is there a recommended approach for recovering the full-rank solution?

Looking forward to your insights!

[rangell]

You can use the function reconstruct_from_sketch:

usbs/utils/common.py

Lines 100 to 116 in f36da24

	def reconstruct_from_sketch(
	Omega: Array,
	P: Array,
	approx_eps: float = 1e-6
	) -> Tuple[Array, Array]:
	n = Omega.shape[0]
	rho = jnp.sqrt(n) * approx_eps * jnp.linalg.norm(P, ord=2)
	P_rho = P + rho * Omega
	B = Omega.T @ P_rho
	B = 0.5 * (B + B.T)
	L = jnp.linalg.cholesky(B)
	E, Rho, _ = jnp.linalg.svd(
	jnp.linalg.lstsq(L, P_rho.T, rcond=None)[0].T,
	full_matrices=False, # this compresses the output to be rank `R`
	)
	Lambda = jnp.clip(Rho ** 2 - rho, 0, jnp.inf)
	return E, Lambda

[ericzzj1989]

Thank you for your response!

I tried using reconstruct_from_sketch, but after applying SVD to the reconstructed matrix, I still only get sketch_dim singular values. Apologies if my previous question wasn’t clear.

What I actually want is not just a low-rank approximation, but rather a way to obtain the full n-dimensional optimal solution while still benefiting from sketching. I understand that turning off sketching would allow me to recover the full solution, but in my case, X has a very large dimension, so disabling sketching isn’t practical.

Is there a way to obtain n-dimensional eigenvectors while still using sketching? I would really appreciate any insights you might have!

Thanks again for your help!

[rangell]

AFAIK, there isn't a way to obtain a higher-rank approximation. The low-rank reconstruction will give you n-dimensional eigenvectors, but you just won't have more than the sketch-dimension worth of non-zero eigenvalues. You will still get an n-by-n matrix, it will just be the best low-rank approximation possible based on the sketch-dimension.

Let me know if you have any further questions.

[ericzzj1989]

Thank you for the clarification!

In my case, suppose X is 31 × 31 and I set sketch_dim = 5. After computing the low-rank reconstruction using:

R, _ = reconstruct_from_sketch(sdp_state_final.Omega, sdp_state_final.P)

I get a matrix R of shape 31 × 5.

I’m wondering how exactly do I obtain the n-dimensional eigenvectors (i.e., the 31-dimensional eigenvectors you mentioned)? Should I take the first left singular vector U[:, 0] from the SVD of R, or is there another recommended approach?

Apologies if this is a basic question, and I really appreciate your patience!

[rangell]

The columns of R (the first value returned from reconstruct_from_sketch) will be your eigenvectors.

reconstruct_from_sketch returns the eigendecomposition of a low-rank approximation of the solution (the first return value is the eigenvectors and the second is the eigenvalues).