That is an excellent answer on so many levels! Thank you for that.

i) Both algorithms (algebraic and matrix form) produces the same output
and the same allocations. Check one.

ii) The matrix QP form runs faster. So matrix form is better to use. Check two.

iii) The [-EV,2*Cov] did not make any sense when I looked at it first but when I saw that the
default form for QP was  c' x + (1/2) x' H x it all made sense. Which means  x' COV x - Ev'x which is exactly
the form we want. Again thank you for your expertise :-)

You only have with(DirectSearch) in the begining which leads to that Search(r, {con1, con2, con3, con4, con5})
produces an ouput of 0.

You have a restart before you simulate the data. If you put with(DirectSearch) after that restart you will get the same error messege as me!

The search procedure does not work for my purpose. It just freezes and when you stop it you get:

with(DirectSearch):
Search(r, {con1, con2, con3, con4, con5});

Warning, initial point [r = .900000000000000022] does not satisfy the inequality constraints; trying to find a feasible initial point
Warning,  computation interrupted

Following Robert solution to only invest a fraction




The expected value is:




We can plot the 3d chart as follows:

@Markiyan Hirnyk Thank you for your suggetion.

I have attached my cone problem MCone.mw

Which procedure in the DirectSearch package v. 2 would you recommend me using to solve the problem most efficiently?

LL_54)_Interior_Poin.mw

@Markiyan Hirnyk. Thank you. I have uploaded the worksheet which contains comments. The comments
where automatically removed when I posted here for some reason.

@Acer you are right. I wanted it to be general and not only solve a specific LP, however the complexity of the task made me give up. In, regards to efficient implimntation. True, there is most likely more efficient code, however I am no expert when it comes to efficient coding so I leave that to some one ells.


The objective was to create an efficient solver that I could use for cone optimization. However, I think I will give up such task because it is to difficult.

When would you need to calculate the "coefficients of a differential polynomial"
Could anyone give a practical example?

The problem "How to count grains of sand in a heap with Maple?" sounds very interesting. However,
I will not look into it because I dont understand anything of what you have done with maple.

If you explain a litlte bit more ie why do you use a system like {x+2*y+5*z-10^2012<=0, x+3*y+4*z-10^2012<=0}  then maybe you can awake my interest.

I did some further simulations:

restart:
with(plots):
with(Statistics):
randomize():

X := proc (z) local p, W1, W2, i;

p := Sample(Bernoulli(.55), 1000):
W1[0] := 100;
W2[0] := 100;

for i to 1000 do

if p[i] = 1 then W1[i] := 1.4*W1[i-1] else W1[i] := .6*W1[i-1] end if;
if p[i] = 1 then W2[i] := 1.1*W2[i-1] else W2[i] := .9*W2[i-1] end if

end do;

if z = w1 then W1[1000] else W2[1000] end if :

end proc:

"Expected Value-1" = ExpectedValue([seq(X(w1), i = 1 .. 1000)]);
"Expected Value-2" = ExpectedValue([seq(X(w2), i = 1 .. 1000)]);

                     "Expected Value-1" = 0.02020255810
                     "Expected Value-2" = 9.081383444 10


It would also be interesting to plot the saddle equilibrium in a 3d-space.
For example we have seen previosusly that our expected value will be largest
the smaller the bet size is EV(1.4, 0.4)> EV(1.1, 0.9)

I can plot it in a 2D-plane as follows:

plot(.55*log((1+c)*100)+.45*log((1-c)*100), c = 0 .. .5)

Any Idea how to plot it in 3D-space with time included?

For the first 3 time periods we get:

restart:
for i to 3 do
W1[i] := .55*(1+c)*W1[i-1]+.45*(1-c)*W1[i-1]
end do ;


W1[1] := 0.55 (1 + c) W1[0] + 0.45 (1 - c) W1[0]
      
W1[2] := 0.55 (1 + c) (0.55 (1 + c) W1[0] + 0.45 (1 - c) W1[0])

          + 0.45 (1 - c) (0.55 (1 + c) W1[0] + 0.45 (1 - c) W1[0])

W1[3] := 0.55 (1 + c) (0.55 (1 + c) (0.55 (1 + c) W1[0] + 0.45 (1 - c) W1[0])

   + 0.45 (1 - c) (0.55 (1 + c) W1[0] + 0.45 (1 - c) W1[0])) + 0.45 (1 - c)

  (0.55 (1 + c) (0.55 (1 + c) W1[0] + 0.45 (1 - c) W1[0])

   + 0.45 (1 - c) (0.55 (1 + c) W1[0] + 0.45 (1 - c) W1[0]))

I would like to see an application with an interior point method that can
solve a second order cone problem really really fast like mosek. I have done some
work on the Affine Scaling Algorithm which can be found here: LL_54)_Interior_Poin.mw
But I am unsure if it will work for a non-linear constraint like a cone?!

Cholesky, excellent stuff !

I also figured out how to do it in another way. Let R be the upper triangular factor of the Cholesky Decomposition matrix. In Maple LUDecomposition(Cov, 'method' = 'Cholesky', output = 'U'). Then we have:




S=R'*R   w'.R'=(R.w)'    sqrt(x'x)=norm(x)

w'.S.w < r^2     w'.R'.R.w< r^2    (R.w)'.(R.w)< r^2
sqrt((R.w)'.(R.w))< sqrt(r^2)    norm(R.w)< r

Cholesky, excellent stuff !

I also figured out how to do it in another way. Let R be the upper triangular factor of the Cholesky Decomposition matrix. In Maple LUDecomposition(Cov, 'method' = 'Cholesky', output = 'U'). Then we have:




S=R'*R   w'.R'=(R.w)'    sqrt(x'x)=norm(x)

w'.S.w < r^2     w'.R'.R.w< r^2    (R.w)'.(R.w)< r^2
sqrt((R.w)'.(R.w))< sqrt(r^2)    norm(R.w)< r

Ok, thank you for that.
Yes maybe R is not equal to R' since the cholesky decomp matrix is trianglular.
The covariance matrix S on the other hand is symmetrical hence
S=S' and S=R*R' where R is the Cholesky decomposition matrix
S'=(R*R')'

but this is as far as I got....how can you rewrite S'=(R*R')' so it fits in

w'*S*w < r^2     w'*R*R'*w< r^2    w'*R'*R*w<r^2     (R*w)'*(R*w)< r^2
sqrt((R*w)'*(R*w))< sqrt(r^2)    norm(R*w)< r

Ok, thank you for that.
Yes maybe R is not equal to R' since the cholesky decomp matrix is trianglular.
The covariance matrix S on the other hand is symmetrical hence
S=S' and S=R*R' where R is the Cholesky decomposition matrix
S'=(R*R')'

but this is as far as I got....how can you rewrite S'=(R*R')' so it fits in

w'*S*w < r^2     w'*R*R'*w< r^2    w'*R'*R*w<r^2     (R*w)'*(R*w)< r^2
sqrt((R*w)'*(R*w))< sqrt(r^2)    norm(R*w)< r

On a separate note but on the same topic

This strategy does not work since you will sooner or later always go bancrupts but it is interesting to do some calculations any how.

You have three dozen bets on the roulette table; 1-12, 13-24, 25-36
If you correctly select the right one you will win 2*bet size.

We can now say that we always want to win 1.
The system we need to solve after 10 loses is;

print(B[1] = Initial*Bet*Size);
for i from 2 to 10 do
print(2*B[i]-add(B[j], j = 1 .. i-1) = 1)
end do;

                          B[1] = Initial Bet Size
                              2 B[2] - B[1] = 1
                          2 B[3] - B[1] - B[2] = 1
                       2 B[4] - B[1] - B[2] - B[3] = 1
                   2 B[5] - B[1] - B[2] - B[3] - B[4] = 1
                2 B[6] - B[1] - B[2] - B[3] - B[4] - B[5] = 1
            2 B[7] - B[1] - B[2] - B[3] - B[4] - B[5] - B[6] = 1
         2 B[8] - B[1] - B[2] - B[3] - B[4] - B[5] - B[6] - B[7] = 1
     2 B[9] - B[1] - B[2] - B[3] - B[4] - B[5] - B[6] - B[7] - B[8] = 1
 2 B[10] - B[1] - B[2] - B[3] - B[4] - B[5] - B[6] - B[7] - B[8] - B[9] = 1

We can solve for the optimal bet sizes as follows:

restart;
z := [x[1] = 10, seq(2*x[i]-add(x[j], j = 1 .. i-1) = 1, i = 2 .. 14)];
evalf(solve(z, {seq(x[i], i = 1 .. 14)}));
assign(%);
-add(x[i], i = 1 .. 14);
-add(x[i], i = 1 .. 13)+2*x[14];

  {x[1] = 10., x[2] = 5.500000000, x[3] = 8.250000000, x[4] = 12.37500000,

    x[5] = 18.56250000, x[6] = 27.84375000, x[7] = 41.76562500,

    x[8] = 62.64843750, x[9] = 93.97265625, x[10] = 140.9589844,

    x[11] = 211.4384766, x[12] = 317.1577148, x[13] = 475.7365723,

    x[14] = 713.6048584}