Poor performance of linprog in practice

I have to solve a dynamic programming problem using a linear programming approach. For details, please see this paper. The LP that I want to solve is:

min c'*v

s.t.

A*v>=u,

where c is n*1, v is n*1, A is n^2*n, u is n^2*1.

The min is with respect to v, the value function of the original DP problem. I have a moderate number of variables, n=300 and m=n^2*n=90000 linear inequalities as constraints. No bound constraints on v.

I use the Matlab function linprog which in turn is based on the solver HIGHS (since R2024a). The code is slow for my purposes (i.e. a brute-force value iteration is much faster). Moreover, linprog gives correct results only if I set the option 'Algorithm','dual-simplex-highs'. With other algorithms, it gets stuck.

After profiling the code, it turns out that the bottleneck is line 377 of linprog:

[x, fval, exitflag, output, lambda] = run(algorithm, problem);

I was wondering if there is a way to speed up the code. Any help or suggestion is greatly appreciated! I put below a MWE to illustrate the problem.

clear,clc,close all
%% Set parameters
crra  = 2;
alpha = 0.36;
beta  = 0.95;
delta = 0.1;
%% Grid for capital
k_ss = ((1-beta*(1-delta))/(alpha*beta))^(1/(alpha-1));
n_k    = 300;
k_grid = linspace(0.1*k_ss,1.5*k_ss,n_k)';
%% Build current return matrix, U(k',k)
cons = k_grid'.^alpha+(1-delta)*k_grid'-k_grid;
U_mat = f_util(cons,crra);
U_mat(cons<=0) = -inf;
%% Using LINEAR PROGRAMMING
% min c'*v
% s.t.
% A*v>=u, where c is n*1, v is n*1, A is n^2*n, u is n^2*1
n     = length(k_grid);
c_vec = ones(n,1);
u_vec = U_mat(:); %% U(k',k), stack columnwise
%% Build A matrix using cell-based method
tic
A = cell(n,1);
bigI = (-beta)*speye(n);
for i=1:n
    temp = bigI;
    temp(:,i) = temp(:,i)+1;
    A{i} = temp;
end
A = vertcat(A{:});
disp('Time to build A matrix with cell method:')

Time to build A matrix with cell method:

toc

Elapsed time is 0.021824 seconds.

%% Call linprog
% 'dual-simplex-highs' (default and by far the best)
options = optimoptions('linprog','Algorithm','dual-simplex-highs');
tic
[V_lin,fval,exitflag,output] = linprog(c_vec,-A,-u_vec,[],[],[],[],options);

Optimal solution found.

disp('Time linear programming:')

Time linear programming:

toc

Elapsed time is 1.373489 seconds.

if exitflag<=0
    warning('linprog did not find a solution')
    fprintf('exitflag = %d \n',exitflag)
end
%% Now that we solved for V, compute policy function
RHS_mat = U_mat+beta*V_lin; % (k',k)
[V1,pol_k_ind] = max(RHS_mat,[],1);
pol_k = k_grid(pol_k_ind);
% Plots
figure
plot(k_grid,V1)

figure
plot(k_grid,k_grid,'--',k_grid,pol_k)

function util = f_util(c,crra)
util = c.^(1-crra)/(1-crra);
end

PROFILE

[Torsten]

What is "brute-force value iteration" ?

[Matt J]

I have to solve a dynamic programming problem using a linear programming approach.

What do you mean by you "have to"? If speed is a priority, it seems like it would be better to solve it with a dynamic programming algorithm.

[Alessandro]

Thanks for your comments!

@Torsten: I have this DP problem to solve:

for all i=1,..,n, where n denotes the dimension of the state space. V is the value function, U is a current-period payoff. Brute-force value iteration means that I solve the Bellman equation above with fixed point iterations on a discrete grid. I guess V0(i) for all i=1,..,n, I plug V0 into the right-hand side of the equation, I maximize over j (this is the very costly step!!) and I get a new value function, let's call it V1. If ||V1-V0||<tolerance, I stop, otherwise I keep on iterating until I find the fixed point.

@Matt J I don't have to in a literal sense. There are some papers in the computational literature that show that if you recast a DP problem as LP problem you can achieve significant speed gains. See e.g. the paper that I have linked in my OP. I am surprised that in this simple example the LP approach is slower than value function iteration. I guess it might depend on the solver that Matlab uses behind the scenes. Trick and Zin (1993) use CPLEX which is different from HIGHS used by Matlab in linprog. They show that LP can be much faster than value iteration (see graph below).

Answer 1

[Matt J]

There are some papers in the computational literature that show that if you recast a DP problem as LP problem you can achieve significant speed gains. ... Trick and Zin (1993) use CPLEX which is different from HIGHS using by Matlab in linprog. They show that LP can be much faster than value iteration (see graph below).

You could try using CPLEX. However, I don't think a 1993 paper is going to be remotely applicable to modern computation. In 1993, there were no GPUs or even multi-core CPUs. You would be very privileged if you had more than 16 MB RAM.

An implemention of DP value iteration with today's computing resources is going to be highly parallelizable with the Parallel Computing Toolbox - much more so than linprog.

[Alessandro]

Thanks for your answer. Maybe I should not have cited the paper from 1993, but nowadays there are parallel implementations of the dual simplex algorithm for LP, see e.g. https://ergo-code.github.io/HiGHS/dev/#Solution-algorithms and https://link.springer.com/article/10.1007/s12532-017-0130-5. If the Matlab linprog is based on HIGHS (as it seems), it should do some parallelization in the background.

The availability of GPUs and multi-core CPUs speeds up DP value iteration but also the LP method, so I don't see why the LP approach would become obsolete for this type of problems. The issue is that linprog does not perform very well

[Matt J]

The availability of GPUs and multi-core CPUs speeds up DP value iteration but also the LP method, so I don't see why the LP approach would become obsolete for this type of problems.

Because value iteration can be written as a sequence of matrix multiplications and other elemental matrix operations, which Matlab is super-well optimized for. Conversely, linprog requires the allocation of memory for all kinds of submatrices and tableaus and doesn't really have a predictable memory access pattern. It really doesn't surprise me that it runs faster than linprog.

However, if you show your "brute force" implementation of value iteration, the community might be able to examine whether you've implemented that in the fastest possible way.

[Matt J]

Also, I notice you don't have a stochastic state transition in your DP problem. That simplifies value iteration greatly, and makes it extra fast. You can see below that one step of value iteration is only 1e-4 seconds. You could do 1000 iterations without breaking a sweat.

U=rand(300);
V=rand(300,1);
beta=0.5;

timeit(@()valueIter(U,V,beta))

ans = 8.0827e-05

function  V=valueIter(U,V,beta)
 V=max(U+beta*V',[],2);
end