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

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)];

[David Hill]

I don't understand your output, it repeats (1*10+3)*12+5; (1*10+3)*12+5; twice in the middle. Is that correct?

[L Smith]

Thanks David. It was an error and I have just updated the questions. I want to calculate C such that:

C = [1;1*10+3;(1*10+3)*12+5; ((1*10+3)*12+5)*14+7] = [1;13;161;2261]

[Stephen23]

Your example output repeats 1*10+3 twice in the middle (elements two and three), so your formula has five terms but your example numeric vector has only four terms. Is that correct?

[Matt J]

There is an inconsistency. This expression,

C=[A(1);A(1)*B(2)+A(2);(A(1)*B(2)+A(2))*B(3)+A(3);((A(1)*B(2)+A(2))*B(3)+A(3))*B(4)+A(4);...
   (((A(1)*B(2)+A(2))*B(3)+A(3))*B(4)+A(4))*B(5)+A(5);((((A(1)*B(2)+A(2))*B(3)+A(3))*B(4)+A(4))*B(5)+A(5))*B(6)+A(6);...
   (((((A(1)*B(2)+A(2))*B(3)+A(3))*B(4)+A(4))*B(5)+A(5))*B(6)+A(6))*B(7)+A(7);...
   ((((((A(1)*B(2)+A(2))*B(3)+A(3))*B(4)+A(4))*B(5)+A(5))*B(6)+A(6))*B(7)+A(7))*B(8)+A(8)];

doesn't use B(1) anywhere, while this example,

   C = [1;1*10+3;1*10+3; (1*10+3)*12+5; ((1*10+3)*12+5)*14+7 ] 

uses B(1) but doesn't use B(end).

[L Smith]

Thank you Matt.

This is a typo as I was copying large number of elements to explain what the loop should do. So it should contain B(1) in the second row and B(end) should not be used. This means I have to correct each successive rows, so I hope I have done it correclty this time.

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)];

Earlier you suggested using one of your functions on Mathworks website to convert the operation generating C to matrix form. In your opinion, between the two methods (i.e. using your function and the code you provided below) is more efficient in terms of run-time or they are pretty much same?

[Matt J]

The version here will be more efficient, but for such a small problem as this, the difference in efficiency will be irrelevant. Also, note that the problem-based optimization framework is slower than the solver-based. You are already working in a framework that does not prioritize speed. But again that's fine, because a problem this small won't benefit from speed optimization.

[L Smith]

Thank you Matt

My actual problem is of much larger scale than the code I shared as I thought it will be easier for anyone to help out if I share a simple script that has all the data with a specific problem.

You are right that problem based approach is much slower than the solver based. I did a simple check on this example and run time is 0.099898 seconds using the Problem based approach, whereas it took only 0.029213 seconds using the Solver-based approach.

I would have preferred coding the problem using the Solver-based approach but I struggle to work out how to code the constraint in the form A.x <= b when I have operations across different elements of matrices. I used prob2struct to run it using the Solver-based approach. Looking at f.Aineq it shows a 8x5 double sparse matrix which I don't know how to code using the Solver-based approach.

Not sure if it is possible to generates MATLAB code using prob2struct to see how this problem will be formulated in the solver-based approach apart from seeing the final result.

[Matt J]

My actual problem is of much larger scale than the code I shared as I thought it will be easier for anyone to help out if I share a simple script that has all the data with a specific problem.

The problem structure you have shared will quickly encounter numerical difficulties as you try to generalize it to larger numbers of variables. Even for the slight increase in dimension shown in the following example data, you can see that the matrix elements blow up very quickly and will ultimately overflow.

A=[1;3;5;7;9;11] ; B = [10;12;14;19;21;23]

d=cumprod(B(1:end-1));
d=[1;d(:)]; 
M=tril(d./d.')

M =

           1           0           0           0           0           0
          10           1           0           0           0           0
         120          12           1           0           0           0
        1680         168          14           1           0           0
       31920        3192         266          19           1           0
      670320       67032        5586         399          21           1

[L Smith]

The matrix B will not have more than 100-120 elements but decision variables which are 5 in the example could be around 10,000 (one dimensional). Is there a way of making it efficient to avoid the numerical problems that you have mentioned?

[Matt J]

I can't see how the problem generalizes when A and B are of different sizes.

[L Smith]

Apologies for not being clear in my question. I don't mean that the dimension of A and B will be different. But number optimization variables 'bonds' as in the code below could be 10,000. So lets say B has 100 rows, A will also have 100 rows. However, in terms of solution, bonds will give 10,000 values.

I think if this results in a 100x100 matrix for M, it should not overflow. Just out of interest if dimension of A and B both was to be much larger than this, is it possible to reformulate it so it does not overflow? I can't think of any way around this.

A=-(cashFlows*bonds-obligations);
d=cumprod(B(1:end-1));
d=[1;d(:)]; 
M=tril(d./d.');

[Matt J]

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

[L Smith]

B(i) in most cases (but not always) be greater than 1. 2^100 seems quite a big value

[Matt J]

That's if the B(i) are all less than 1.1, asin you rexample, it should not be a big issue.

Answer 1

[Matt J]

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)

[Matt J]

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)

[L Smith]

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.

Answer 2

[Stephen23]

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

[L Smith]

Thanks Stephen for the code

You are right that I repeated a term "1*10+3" that was an error. I have corrected it now:

C = [1;1*10+3;(1*10+3)*12+5; ((1*10+3)*12+5)*14+7] = [1;13;161;2261]

I think it may be possible to vectorise it but it is likely to be difficult. I believe it will probably involve multiplication by some upper triangle matrices of 1 with some other operations but I cannot work it out.

The reason I want to avoid a loop is simply becuase I will be using it in an optimization problem (where one of the matrices A or B will have the optimzation variables) and loops can be quite slow.

I will wait for others to respond before accepting your answer.

[Stephen23]

"The reason I want to avoid a loop is simply becuase ... loops can be quite slow."

This is a classic MATLAB myth.

Loops in themselves are not slow, and for many versions of MATLAB the JIT engine has optimized their internal operations. A well designed loop is fast, a badly designed loop is slow: this has nothing to do with the loop per se, but the code design. For that matter it is also possible to write very inefficient vectorized code (just because it is vectorized does not mean that somehow it is magically efficient).

Note that vectorized code can be very slow, especially when it requires large intermediate arrays to be stored in memory (and which would proably be required for this task).

Also keep in mind that the vectorized code will probably not be simple (as your own attempts have no doubt shown you), thus your code's efficiency will be:

runtime           % unknown, possibly slower
writing time      % slower
debugging time    % slower
maintenance time  % slower

[L Smith]

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.

[Daniel M]

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.

[L Smith]

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