Sum of Inverse of Variables in an Optimization Problem

Question

I have the following optimization problem: $$ \begin{array}{ll} \text{Minimize} & \frac{1}{x_1} + \frac{1}{x_2} + \ldots + \frac{1}{d_n} \\ \text{Subject to} & A x \leq b \end{array} $$ where $x = [x_1~~x_2~~\ldots~~x_n]^T$, $x_1,x_2,\ldots,x_n \in \left\{1,\frac{1}{2},\frac{1}{3},\ldots, \frac{1}{k} \right\}$for some $k$, and $A$ and $b$ have positive entries ($A$ is big but actually very sparse).

The dimension $n$ is not too large (I'd say at most 100), but the number of constraints can be huge (hundreds of thousands). They're constraints relating pairs of $x$'s. How can one attempt solving such a problem? Is there a way of transforming this into some standard form that can be solved using some optimization package like Mosek?

Lastly, I can get away with solving a continuous relaxation of this problem if the discrete version proved to be very hard. If the variables are continuous, the feasible space becomes a standard convex polytope, but the objective function remains problematic. Any ideas are appreciated.

Thank you.

Answer 1

For the discrete version, it can be cast as a mixed-integer linear program. You just have to note that every element $x_i$ can be written as $x_i = \sum_{j=1}^k \frac{\delta_{ij}}{j}$ where $\sum_{j=1}^k \delta_{ij} = 1$ where $\delta_{ij}$ is binary. Using the same binary variables the inverses are simply $x_{i}^{-1} = \sum_{j=1}^k j\delta_{ij}$

Here is the YALMIP code to setup and solve such a problem as a proof-of-concept (disclaimer, the MATLAB Toolbox YALMIP, which interfaces Mosek, is developed by me)

k = 5;
n = 10;
delta = binvar(n,k,'full')
x =    delta*(1./(1:k))';
xinv = delta*(1:k)';
A = randn(2*n,n);
b = 10*rand(2*n,1);
Constraints = [A*x <= b, sum(delta,2) == 1];
Objective = sum(xinv);
optimize(Constraints,Objective)
value(x)

The continuous version you speak of is a convex program, as the inverse is convex on the positive orthant. An inverse scalar can be reformulated to a model involving second order cones rather easily, and thus you can solve this as a second-order cone program in Mosek. YALMIP does that automatically if you use the convexity-aware cpower operator.

x = sdpvar(n,1);    
Constraints = [A*x <= b, 0 <= x <= 1];
optimize(Constraints,sum(cpower(x,-1)))
value(x)

Manually, you introduce new variables $t_i$ to replace the inverses and add the constraints $x_i^{-1} \leq t_i$, i.e., $1 \leq x_it_i$, which can be written as the cone constraint $\Big|\Big| \begin{matrix}1\\t_i-x_i\end{matrix}\Big|\Big|\leq x_i+t_i$

Now, whether this scales for large instances, that's a completely different question.