I have a general square matrix [[1,1,0],[1,1,0],[0,0,0]] (the element is on a binary field--Galois Fields) I want to find the rank of this matrix (like the function: np.linalg.matrix_rank). How can I describe this process using gurobi?

Due to the elements in the matrix being related to the variables in gurobi model, I can't use the function of Numpy directly.

I tried asking a question in the Gurobi community before and got a small model that worked. This model has an objective function to get the rank.

When I embedded it in my original big model (which already has an objective function), the original model affected the optimization process of the rank-finding process. That is, in this multi-objective optimization problem, the rank is always wrong.

This would render the entire model unusable since the rank should be a deterministic property of the matrix. But now because the optimization process can't get the best, the rank property becomes uncertain.

My original goal was to maximize a "GObj" that is >= 0. If this value equals 1, it means the model has a solution. In the big model, the rank is required several times (several rank-finding processes are also related to each other).

For the rank process to be optimal, I tried to adjust the priority of the object, but even though I set the priority of the rank-finding very high, I still got a lot of wrong ranks. Similarly, I try not to optimize the "GObj", just limit it to >=1 so that the model can have a solution. I still got a lot of wrong ranks. I also tried to make a lot of small adjustments, but nothing worked.

I think an easy-to-understand solution is to figure out how to describe the process of finding the rank of matrix without objective function. So how can I describe it without objective function?

If I don't change the way the rank is found, how can the modeling be done so that the rank-finding procedure is optimal while maximizing "GObj"?

To describe my problem in detail, I wrote a demo and uploaded it to Github: https://github.com/wenny-kiki/test_math


1 Answer 1


In theory we can use the fact that rank = the number of non-zero singular values in a singular value decomposition (SVD).

import gurobipy as gp
import numpy as np

m = gp.Model()

p,q = 4,3 # as an example, overwritten for testing

# The matrix we want to calculate rank of
M = m.addMVar((p,q), lb=-float("inf"), name="M")

# # uncomment the following lines for testing
# A = np.random.randint(-2,2, size=(p,q))
# m.addConstr(M == A)

U = m.addMVar((p,p), lb=-float("inf"), name="U")

ub = np.zeros((p,q))
np.fill_diagonal(ub, np.inf)
lb = np.zeros((p,q))
np.fill_diagonal(lb, -np.inf)
S = m.addMVar((p,q), lb=lb, ub=ub, name="S")

V = m.addMVar((q,q), lb=-float("inf"), name="V")

m.addConstr([email protected] == np.eye(p))
m.addConstr([email protected] == np.eye(q))

aux = m.addMVar((p,q), lb=-float("inf"), name="aux")
m.addConstr(U@S == aux)
m.addConstr(aux@V == M)

nzpos = m.addVars(min(p,q), vtype="B")
nzneg = m.addVars(min(p,q), vtype="B")
for i in range(min(p,q)):
    m.addConstr((nzpos[i]==1)>>(S[i,i] >= 0.0001))
    m.addConstr((nzpos[i]==0)>>(S[i,i] <= 0.00009))
    m.addConstr((nzneg[i]==1)>>(S[i,i] <= -0.0001))
    m.addConstr((nzneg[i]==0)>>(S[i,i] >= -0.00009))

rank = m.addVar()
m.addConstr(rank == nzpos.sum() + nzneg.sum())

In practice this is not going to work though, and Gurobi will complain that "Model contains variables with very large bounds participating in product terms". We can introduce bounds on U (and V) but in doing so requires $UU^T = I$ to be replaced with $UU^T = aI$ for some non-zero scalar $a$ (and similarly for V). What you get is a model which will probably only be solvable for tiny matrices - not to mention be prone to numerical issues and hence incorrect results.

  • $\begingroup$ This is a nice direction to try.I tried to run the code with a 3 x 3 square, and it took too long to run. I can't even get the run result to test if it's correct. $\endgroup$
    – wenny
    Apr 7 at 7:42
  • $\begingroup$ I think fortunately this is actually possible in theory.I wonder if restricting the process to the binary field would make the model easier?I set the variables as follows,but the running result shows that there is no solution,I don't know what's wrong with that?Thanks for your answer! ``` M = m.addMVar((p, q), vtype="B", name="M") U = m.addMVar((p, p), vtype="B", name="U") modulo_U = m.addMVar((p, p), lb=-5, ub=5, vtype="I", name="modulo_U") m.addConstr(U @ U.T == np.eye(p) + 2*modulo_U) m.addConstr(V @ V.T == np.eye(q) + 2*modulo_V) ``` $\endgroup$
    – wenny
    Apr 7 at 8:03
  • $\begingroup$ Not sure whether the SVD result translates to the binary field or not. Maybe ask in math.stackexchange? $\endgroup$
    – Riley
    Apr 7 at 11:47
  • $\begingroup$ Thanks. I'll try. I searched the web and asked chatgpt before trying to modify the code you gave. chatgpt says SVD works on binary field, but I'm really not sure yet. $\endgroup$
    – wenny
    Apr 8 at 5:26
  • $\begingroup$ Also, I had tried to modify the code to include bounds on the matrices which are multiplied together, which requires a few auxiliary MVars and constraints to be introduced as I alluded to, and it did help. I seemed to be able to get the rank for 3x3 matrices without too much trouble, but sizes like 5x5 seemed to struggle. $\endgroup$
    – Riley
    Apr 8 at 12:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.