MATLAB Answers

L Smith
1

Optimization with loops Error in optim.prob​lemdef.Opt​imizationP​roblem/sol​ve

Asked by L Smith
on 7 Oct 2019
Latest activity Edited by Matt J
on 11 Oct 2019
Q1. I want to perform certain operation on two large column vectors and calculate a cumulative sum in a certain order without using loops. For simplicity assume that I have two column vectors:
A=[1;3;5;7] and B = [10;12;14;19]
and I want to find a vector C = [1;1*10+3;1*10+3; (1*10+3)*12+5; ((1*10+3)*12+5)*14+7 ] = [1;13;161;2261]
What will be the MATLAB code to calculate vector C in MATLAB without using loops?
Q2. I have provided my code below and I have used the helpful suggestion given in the answer to use a loop in an optimization problem. The use of a optimization expression C in the code below is intended to be similar in terms of relevant elements to add and multiply as in the example given above but it contains the optimization variables. I get the error messages shown below when I run the code. I have also tried to write operations on C manually as shown in the last line of the code but I get the same error messages at the last line of the code Solution = solve(ALM).
How can I correct the code?
ERROR MESSAGES
Index exceeds the number of array elements (0).
Error in optim.internal.problemdef.SubsasgnExpressionImpl/computeLinearCoefficients
Error in optim.internal.problemdef.ExpressionTree/extractLinearCoefficients
Error in optim.internal.problemdef.ExpressionForest/extractLinearCoefficients
Error in optim.problemdef.OptimizationExpression/extractLinearCoefficients
Error in optim.problemdef.OptimizationConstraint/extractLinearCoefficients
Error in optim.problemdef.OptimizationProblem/compileConstraints
Error in prob2structImpl
Error in optim.problemdef.OptimizationProblem/solve
prices = [99.74 91.22 98.71 103.75 97.15];
cashFlows = [4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 2.5 5 4; 4 5 102.5 5 4;4 5 0 105 104;4 105 0 0 0; 104 0 0 0 0];
obligations = [5 7 7 6 8 7 20 0]'*1000;
nt=size(cashFlows,1)
nb=size(cashFlows,2)
Rates = [0.01; 0.015; 0.017;0.019;0.02;0.025;0.027;0.029];
EndTimes = (1:nt)';
Disc = rate2disc(1,Rates,EndTimes);
%Number of bonds available
nBonds = [10;100;20;30;5]
ALM = optimproblem;
bonds = optimvar('bonds',nb,'Type','integer','LowerBound',0,'UpperBound',nBonds);
ALM.ObjectiveSense = 'minimize';
ALM.Objective = prices*bonds;
%Define the constraint
B=[1.01;1.020024752;1.02101183;1.02502363;1.024009823;1.050370059;1.039082218;1.043109481];
A=-(cashFlows*bonds-obligations);
C = optimexpr(nt,1);
C(1) = A(1);
for k = 2:numel(A)
C(k) = C(k-1) * B(k) + A(k);
end
ALM.Constraints.Const1 = (C.*Disc)/53.844 <=0.05;
%Solve the problem
Solution = solve(ALM)
%Here I show the intended operation on Optimization Expression C (in case
%for loop above does not operate as intended but I get the same error messages
C=[A(1);A(1)*B(1)+A(2);(A(1)*B(1)+A(2))*B(2)+A(3);(A(1)*B(1)+A(2))*B(2)+A(3))*B(3)+A(4);...
((A(1)*B(1)+A(2))*B(2)+A(3))*B(3)+A(4))*B(4)+A(5);
(((A(1)*B(1)+A(2))*B(2)+A(3))*B(3)+A(4))*B(4)+A(5))*B(5)+A(6);...
((((A(1)*B(1)+A(2))*B(2)+A(3))*B(3)+A(4))*B(4)+A(5))*B(5)+A(6))*B(6)+A(7);...
(((((A(1)*B(1)+A(2))*B(2)+A(3))*B(3)+A(4))*B(4)+A(5))*B(5)+A(6))*B(6)+A(7))*B(7)+A(8)];

  14 Comments

I think if this results in a 100x100 matrix for M, it should not overflow.
Are B(i) always integers greater than 1? If so, then M(100,1)=prod(B) will always be at least
>> 2^100
ans =
1.2677e+30
B(i) in most cases (but not always) be greater than 1. 2^100 seems quite a big value
That's if the B(i) are all less than 1.1, asin you rexample, it should not be a big issue.

Sign in to comment.

2 Answers

Answer by Matt J
on 10 Oct 2019
Edited by Matt J
on 10 Oct 2019
 Accepted Answer

I believe it will probably involve multiplication by some upper triangle matrices of 1 with some other operations but I cannot work it out.
Here it is,
d=cumprod(B(1:end-1));
d=[1;d(:)];
M=tril(d./d.');
ALM.Constraints.Const1 = ((M*A).*Disc)/53.844 <=0.05;
Solution = solve(ALM)

  2 Comments

I don't know why the solver has trouble when a for-loop is used, but it isn't strictly related to the for-loop. For example, this version runs without error messages,
B=[1.01;1.020024752;1.02101183;1.02502363;1.024009823;1.050370059;1.039082218;1.043109481];
A=-(cashFlows*bonds-obligations);
C = optimexpr(nt,1);
C(1) = A(1);
for k = 2:numel(A)
C(k) = C(k-1) * B(k) + A(k);
end
ALM.Constraints.Const1 = C <= 5e6;
Solution = solve(ALM)
Thank you Matt.
I think this should be posted as a separate question as it appears to me an issue with the problem-based approach. Another reason why solver-based approach is better.

Sign in to comment.


Answer by Stephen Cobeldick on 8 Oct 2019
Edited by Stephen Cobeldick on 9 Oct 2019

I don't see a simple vectorized solution, but one for loop will be quite efficient:
A = [1;3;5;7]
B = [10;12;14;19]
C = nan(size(A)); % preallocate
C(1) = A(1);
for k = 2:numel(A)
C(k) = C(k-1) * B(k-1) + A(k);
end
Giving:
C =
1
13
161
2261

  5 Comments

Thank you Stephen. You are right that poorly vectorized code can be slow than a loop. I am not a regular user of MATLAB, hence the question.
I have provided my code now in the question. I have not been to make the code work using the code you kindly provided. I believe MATLAB allows only limited set of operations on optimization expressions.
I have a hunch there is a very clever solution using something like conv2 or filter, but it would be much more complicated to understand than the loop.
MATLAB only allows a limited set of operations on optimization expressions, and conv2 and filter are not listed in the documentation here: https://uk.mathworks.com/help/optim/ug/supported-operations-on-optimization-variables-expressions.html

Sign in to comment.