Add absolute value and power operator for the BinaryQuadraticModel

[angystallone]

Currently, only multiplication, addition, subtraction and division of BinaryQuadraticModel by a number are allowed. However, to build a more complex objective function, additional operations are needed. Among them, the absolute value is essential for the definition of the objective as the L1-norm, as in the following example:

# Add binary variables
x = [dimod.Binary(i) for i in range(N)]

# Loop over rows of matrix A
for i in range(rows):
    # Initialize expra
    expra = 0*len(x)
    for j in range(N2):
        if A[i] != 0:
            # Build the linear expression
            expra += A[i]*x[j]
            # Take the abs value
            expraR[i] = abs(expra - number)

# Construct the CQM
cqm = CQM()

## Add the objective (L1-norm)
cqm.set_objective(sum(expraR))

In the code above, for each row of a matrix A, a linear mathematical expression (expra in the script) is formed as a function of the binary variable x. There are no issues in generating expra, but there is a problem when defining expraR, since I cannot calculate the absolute value of the difference inside the parenthesis (L1-norm).

Using this alternative solution does not work:
expraR[i] = ((expra - number) ** 2) ** 0.5
as the power operator is not supported either, returning the following error:
TypeError: unsupported operand type(s) for ** or pow(): 'BinaryQuadraticModel' and 'float'

[arcondello]

Hi @angystallone,

I don't have the full context of your code, but it looks like you're trying to minimize the distance between a bqm and single number. In that case, generally the best way to do it would be something like objective = (bqm - number)**2. You don't actually need to take the square root at all because the CQM solver handles quadratic expressions in addition to linear.

For the feature request, it is often best to think of bqms as polynomials rather than as matrices. If you do that, you'll see why we have not implemented these methods for binary quadratic models in general.

Absolute Value

Taking a simple linear bqm, $E(x, y) = x - y$, let's say that we want to calculate the absolute value of it. $abs(E(x, y))$. We could fairly easily implement an absolute value method on the BQM, probably defined something like

def __abs__(self) -> BinaryQuadraticModel:
    for v, bias in self.iter_linear():
        self.set_linear(v, abs(bias)
    for u, v, bias in self.iter_quadratic():
        self.set_quadratic(u, v, abs(bias))
    self.offset = abs(self.offset)

Let's call the output of this new function $E'(x, y)$ and see that $E(x, y) = x + y$. However, this function would not actually give us what we want. Consider the value of the equation for all possible assignments of $x$ and $y$.

x, y	E(x, y)	abs(E(x, y))	E'(x, y)
0, 0	0	0	0
0, 1	-1	1	1
1, 0	1	1	1
1, 1	0	0	2

Square Root

Similarly, we could implement a square root function. Let's try to take the square root of $(x - y)^2 = x^2 - 2xy + y^2$.
In this case the "square root" is not unique! It could be either be $x-y$ or $y-x$.

The single variable case
In your case it seems that you want to take the absolute value or the square root of binary quadratic models with only a single variable. We could implement something like

def __abs__(self) -> BinaryQuadraticModel:
    if self.num_variables > 1:
        raise TypeError("cannot take the absolute value of bqms with more than one variable")
    for v, bias in self.iter_linear():
        self.set_linear(v, abs(bias)
    self.offset = abs(self.offset)

but that seems like a very specific edge case, and to the point that it would almost be misleading

[angystallone]

Thanks for the quick and clear answer @arcondello.

I see your point about matrices and polynomials. Indeed, expra needs to be defined as a polynomial: the goal here is to first build a linear expression as a function of the variable x for each row of matrix A, and then to minimize the total misfit defined as the sum of these linear expressions relative to a number.

The absolute value of (expra - number) is needed since the objective function must be defined as the L1-norm, and not the L2-norm. Reason for that is I need to map a problem already solved in the classical way onto the D-Wave solver. Unfortunately, in order to reproduce the results, I am not free to make any changes to the problem statement.

More on the problem context: it falls into the category of knapsack problems and in the classical approach the linear expression has been coded by the “model” python module in docplex.mp ([http://ibmdecisionoptimization.github.io/docplex-doc/mp/docplex.mp.model.html]).

[angystallone]

Hi @arcondello, I have an update from my side.

Even switching to L2-norm (so substituting expraR[i] = abs(expra - number) with expraR[i] = (expra - number) ** 2) does not work. This time, I do not get any error message, but the code is stuck, running forever.

[arcondello]

Hi @angystallone , would you be able to share your code? Or, even better, a minimal reproducible example?

[angystallone]

@arcondello Here are the minimal examples for reproducing both the abs value (for L1-norm) and the square (for L2-norm) issues:

from dimod import Binary, CQM
from dwave.system import DWaveSampler, LeapHybridCQMSampler

solver = LeapHybridCQMSampler()
print(solver.properties.keys())

###################################

# Add binary variables
x = [Binary(i) for i in range(100)]

expraR = [[0]*len(x)] * 9
# print (len(expraR),len(expraR[0]))
for i in range(9):
    expra = 0 * len(x)
    for j in range(100):
        expra += 7.23 * x[j]

    ### Taking abs value of the difference is not working, as it is not among the operations that can be applied to a BQM
    # expraR[i] = abs(expra - 470.0)

    # Alternative way to take the abs value, not working too
    expraR[i] = ((expra - 470.0) ** 2) ** 0.5

In this second case, the problem could be due to the size of x (for my code, it proceeds very slowly, then it stops at i = 32):

from dimod import Binary, CQM
from dwave.system import DWaveSampler, LeapHybridCQMSampler

solver = LeapHybridCQMSampler()
print(solver.properties.keys())

###################################

# Add binary variables
x = [Binary(i) for i in range(50000)]

expraR = [[0] * len(x)] * 90
# print (len(expraR),len(expraR[0]))
for i in range(90):
    expra = 0 * len(x)
    for j in range(50000):
        expra += 7.23 * x[j]
    expraR[i] = (expra - 470.0) ** 2
    print(i)

Thanks a lot!