@Stephen Forrest 

i have tested with ([a,b,a,b,a,b], 3)

and successfully to get full combination

this is why i concatenate list 3 times

if this is not the correct way to get full combination, what is the correct way to get full combination and without out of memory?

actually my goal is to group system of polynomials by hilbert series

save hilbert series result in dictionary key, and then save its system of polynomials into its value.

this case is the beginning because it only use 0 and 1 in matrix

more case i want to find is variables's power can be 2 or more. for example

x^2*y .... etc. this will be more than 3 billions.

first thing is to do variables's power can only be 1.

@Mac Dude 

i want to programmatically save the graph of this command

dualaxisplot(listplot(sols[2], color = red), listplot(sols[1], color = blue), style = line, gridlines = false);

save as a a series of png or jpeg by updating input sol[2] and sol[1]

and 

also want to animate the graph by updating sol[2] and sol[1]

@acer 

i use secondary school method to calculate and validate it, i obtain another eigenvector solution

but it is not the same as maple's one.

how to solve by hand in order to be the same as maple?

when doing quantum logic with fuzzy operation and if values is not the same, the result will be different,

if there are many solutions, how do operator theory, von neuman logic consistent among thesis?

InputMatrix3 := Matrix([[31.25,30.8,30.5],[30.8,30.5,0],[30.5,0,0]]);
NewInput3 := MatrixMatrixMultiply(Transpose(InputMatrix3), InputMatrix3);
Old_Asso_eigenvector := Eigenvectors(NewInput3);
FirstEigenValue := solve(Determinant(NewInput3-Matrix([[lambda1, 0, 0], [0, lambda1, 0], [0, 0, lambda1]])), lambda1)[3]; # find back eigenvalue from eigenvector
SecondEigenValue := solve(Determinant(NewInput3-Matrix([[lambda1, 0, 0], [0, lambda1, 0], [0, 0, lambda1]])), lambda1)[2]; # find back eigenvalue from eigenvector
ThirdEigenValue := solve(Determinant(NewInput3-Matrix([[lambda1, 0, 0], [0, lambda1, 0], [0, 0, lambda1]])), lambda1)[1]; # find back eigenvalue from eigenvector
sys1:=NewInput3-Matrix([[FirstEigenValue, 0, 0], [0, FirstEigenValue, 0], [0, 0, FirstEigenValue]]);
sys2:=NewInput3-Matrix([[SecondEigenValue, 0, 0], [0, SecondEigenValue, 0], [0, 0, SecondEigenValue]]);
sys3:=NewInput3-Matrix([[ThirdEigenValue, 0, 0], [0, ThirdEigenValue, 0], [0, 0, ThirdEigenValue]]);
sys1b:=MatrixMatrixMultiply(NewInput3-Matrix([[FirstEigenValue, 0, 0], [0, FirstEigenValue, 0], [0, 0, FirstEigenValue]]),Matrix([[x],[y],[z]]));
sys2b:=MatrixMatrixMultiply(NewInput3-Matrix([[SecondEigenValue, 0, 0], [0, SecondEigenValue, 0], [0, 0, SecondEigenValue]]),Matrix([[x],[y],[z]]));
sys3b:=MatrixMatrixMultiply(NewInput3-Matrix([[ThirdEigenValue, 0, 0], [0, ThirdEigenValue, 0], [0, 0, ThirdEigenValue]]),Matrix([[x],[y],[z]]));
if Rank(sys1) < 3 then
 print("infinite solutions")
else
 print("unique solutions")
end if:
if Rank(sys2) < 3 then
 print("infinite solutions")
else
 print("unique solutions")
end if:
if Rank(sys3) < 3 then
 print("infinite solutions")
else
 print("unique solutions")
end if:
LinearSolve(sys1);
LinearSolve(sys2);
LinearSolve(sys3);
sys1c:=MatrixMatrixMultiply(GaussianElimination(sys1),Matrix([[x],[y],[z]]));
sys2c:=MatrixMatrixMultiply(GaussianElimination(sys2),Matrix([[x],[y],[z]]));
sys3c:=MatrixMatrixMultiply(GaussianElimination(sys3),Matrix([[x],[y],[z]]));
sol1 := solve([sys1b[1,1],sys1b[2,1], z=t], [x,y,z]);
sol1 := solve([sys1c[1,1],sys1c[2,1], z=t], [x,y,z]);
sol2 := solve([sys2b[1,1],sys2b[2,1], z=t], [x,y,z]);
sol2 := solve([sys2c[1,1],sys2c[2,1], z=t], [x,y,z]);
sol3 := solve([sys3b[1,1],sys3b[2,1], z=t], [x,y,z]);
sol3 := solve([sys3c[1,1],sys3c[2,1], z=t], [x,y,z]);
K := Matrix([[subs(t = 1, rhs(sol1[1, 1])), subs(t = 1, rhs(sol1[1, 2])), subs(t = 1, rhs(sol1[1, 3]))],
[subs(t = 1, rhs(sol2[1, 1])), subs(t = 1, rhs(sol2[1, 2])), subs(t = 1, rhs(sol2[1, 3]))],
[subs(t = 1, rhs(sol3[1, 1])), subs(t = 1, rhs(sol3[1, 2])), subs(t = 1, rhs(sol3[1, 3]))]]):
v := ColumnSpace(K);
w1 := v[1][1];
w2 := v[1][2]-(v[1][2].w[1])*w[1]/(w[1].w[1]);
w3 := v[1][3]-(v[1][3].w[1])*w[1]/(w[1].w[1])-(v[1][3].w[2])*w[2]/(w[2].w[2]);
# validate
Validatesol1:=MatrixMatrixMultiply(NewInput3-Matrix([[FirstEigenValue, 0, 0], [0, FirstEigenValue, 0], [0, 0, FirstEigenValue]]),Matrix([[subs(t = 1, rhs(sol1[1, 1]))],[subs(t = 1, rhs(sol1[1, 2]))],[subs(t = 1, rhs(sol1[1, 3]))]]));
Validatesol2:=MatrixMatrixMultiply(NewInput3-Matrix([[SecondEigenValue, 0, 0], [0, SecondEigenValue, 0], [0, 0, SecondEigenValue]]),Matrix([[subs(t = 1, rhs(sol2[1, 1]))],[subs(t = 1, rhs(sol2[1, 2]))],[subs(t = 1, rhs(sol2[1, 3]))]]));
Validatesol3:=MatrixMatrixMultiply(NewInput3-Matrix([[ThirdEigenValue, 0, 0], [0, ThirdEigenValue, 0], [0, 0, ThirdEigenValue]]),Matrix([[subs(t = 1, rhs(sol3[1, 1]))],[subs(t = 1, rhs(sol3[1, 2]))],[subs(t = 1, rhs(sol3[1, 3]))]]));

@Preben Alsholm 

i use solve and validate my own solution, this solution is also correct!!!

but this solution is not the same as maple's one.

how do maple solve for eigenvector?

InputMatrix3 := Matrix([[31.25,30.8,30.5],[30.8,30.5,0],[30.5,0,0]]);
NewInput3 := MatrixMatrixMultiply(Transpose(InputMatrix3), InputMatrix3);
Old_Asso_eigenvector := Eigenvectors(NewInput3);
FirstEigenValue := solve(Determinant(NewInput3-Matrix([[lambda1, 0, 0], [0, lambda1, 0], [0, 0, lambda1]])), lambda1)[3]; # find back eigenvalue from eigenvector
SecondEigenValue := solve(Determinant(NewInput3-Matrix([[lambda1, 0, 0], [0, lambda1, 0], [0, 0, lambda1]])), lambda1)[2]; # find back eigenvalue from eigenvector
ThirdEigenValue := solve(Determinant(NewInput3-Matrix([[lambda1, 0, 0], [0, lambda1, 0], [0, 0, lambda1]])), lambda1)[1]; # find back eigenvalue from eigenvector
sys1:=NewInput3-Matrix([[FirstEigenValue, 0, 0], [0, FirstEigenValue, 0], [0, 0, FirstEigenValue]]);
sys2:=NewInput3-Matrix([[SecondEigenValue, 0, 0], [0, SecondEigenValue, 0], [0, 0, SecondEigenValue]]);
sys3:=NewInput3-Matrix([[ThirdEigenValue, 0, 0], [0, ThirdEigenValue, 0], [0, 0, ThirdEigenValue]]);
sys1b:=MatrixMatrixMultiply(NewInput3-Matrix([[FirstEigenValue, 0, 0], [0, FirstEigenValue, 0], [0, 0, FirstEigenValue]]),Matrix([[x],[y],[z]]));
sys2b:=MatrixMatrixMultiply(NewInput3-Matrix([[SecondEigenValue, 0, 0], [0, SecondEigenValue, 0], [0, 0, SecondEigenValue]]),Matrix([[x],[y],[z]]));
sys3b:=MatrixMatrixMultiply(NewInput3-Matrix([[ThirdEigenValue, 0, 0], [0, ThirdEigenValue, 0], [0, 0, ThirdEigenValue]]),Matrix([[x],[y],[z]]));
if Rank(sys1) < 3 then
 print("infinite solutions")
else
 print("unique solutions")
end if:
if Rank(sys2) < 3 then
 print("infinite solutions")
else
 print("unique solutions")
end if:
if Rank(sys3) < 3 then
 print("infinite solutions")
else
 print("unique solutions")
end if:
LinearSolve(sys1);
LinearSolve(sys2);
LinearSolve(sys3);
sys1c:=MatrixMatrixMultiply(GaussianElimination(sys1),Matrix([[x],[y],[z]]));
sys2c:=MatrixMatrixMultiply(GaussianElimination(sys2),Matrix([[x],[y],[z]]));
sys3c:=MatrixMatrixMultiply(GaussianElimination(sys3),Matrix([[x],[y],[z]]));
sol1 := solve([sys1b[1,1],sys1b[2,1], z=t], [x,y,z]);
sol1 := solve([sys1c[1,1],sys1c[2,1], z=t], [x,y,z]);
sol2 := solve([sys2b[1,1],sys2b[2,1], z=t], [x,y,z]);
sol2 := solve([sys2c[1,1],sys2c[2,1], z=t], [x,y,z]);
sol3 := solve([sys3b[1,1],sys3b[2,1], z=t], [x,y,z]);
sol3 := solve([sys3c[1,1],sys3c[2,1], z=t], [x,y,z]);
K := Matrix([[subs(t = 1, rhs(sol1[1, 1])), subs(t = 1, rhs(sol1[1, 2])), subs(t = 1, rhs(sol1[1, 3]))],
[subs(t = 1, rhs(sol2[1, 1])), subs(t = 1, rhs(sol2[1, 2])), subs(t = 1, rhs(sol2[1, 3]))],
[subs(t = 1, rhs(sol3[1, 1])), subs(t = 1, rhs(sol3[1, 2])), subs(t = 1, rhs(sol3[1, 3]))]]):
v := ColumnSpace(K);
w1 := v[1][1];
w2 := v[1][2]-(v[1][2].w[1])*w[1]/(w[1].w[1]);
w3 := v[1][3]-(v[1][3].w[1])*w[1]/(w[1].w[1])-(v[1][3].w[2])*w[2]/(w[2].w[2]);
# validate
Validatesol1:=MatrixMatrixMultiply(NewInput3-Matrix([[FirstEigenValue, 0, 0], [0, FirstEigenValue, 0], [0, 0, FirstEigenValue]]),Matrix([[subs(t = 1, rhs(sol1[1, 1]))],[subs(t = 1, rhs(sol1[1, 2]))],[subs(t = 1, rhs(sol1[1, 3]))]]));
Validatesol2:=MatrixMatrixMultiply(NewInput3-Matrix([[SecondEigenValue, 0, 0], [0, SecondEigenValue, 0], [0, 0, SecondEigenValue]]),Matrix([[subs(t = 1, rhs(sol2[1, 1]))],[subs(t = 1, rhs(sol2[1, 2]))],[subs(t = 1, rhs(sol2[1, 3]))]]));
Validatesol3:=MatrixMatrixMultiply(NewInput3-Matrix([[ThirdEigenValue, 0, 0], [0, ThirdEigenValue, 0], [0, 0, ThirdEigenValue]]),Matrix([[subs(t = 1, rhs(sol3[1, 1]))],[subs(t = 1, rhs(sol3[1, 2]))],[subs(t = 1, rhs(sol3[1, 3]))]]));

@acer 

actually this reply can be another question, but it has been deleted, so i posted here

below is the steepest descent method to replace the minimze function to prepare write in other programming languages.

expect it has the same result as minimize function in maple

when it evalf a and b, it return 0 only, i do not know whether caused by max and piecewise function this complicated result

with(Optimization):

f1 := -2*x1-x2; f2 := -x1-4*x2; g1 := 2*x1+3*x2-6; g2 := -x1; g3 := -x2;

penalty := lambda1*max(f1-M,0) + lambda2*max(f2-M,0) + (M^2)*(max(g1,0) + max(g2,0) + max(g3,0)):

k := 1:

s := 1:

# write into procedure

L1 := 0.5;

L2 := 0.5;

with(VectorCalculus):

with(LinearAlgebra):

obj := eval(penalty,[lambda1=3,lambda2=0.645,M=1]);

with(Student:-VectorCalculus):

Hf :=Hessian(obj,[x1,x2]);

Hf :=Jacobian(Gradient(obj,[x1,x2]),[x1,x2]);

#Hf := Matrix([[diff(diff(obj, x1), x1), diff(diff(obj, x1), x2)],

#[diff(diff(obj, x2), x1), diff(diff(obj, x2), x2)]]);

G := Matrix([[diff(obj, x1)], [diff(obj, x2)]]);

z1 := ZeroMatrix(1, 2, compact = false);

e := 1;

slow_constant := 8.92;

z1[1, 1] := 0.5;

z1[1, 2] := 0.5;

maxnumber := 100000;

count := 1;

for count from 1 to maxnumber do

if e > 0.001 then

                a := MatrixMatrixMultiply(Transpose(evalf(subs(x2=z1[1, 2],subs(x1=z1[1, 1],G)))), evalf(subs(x2=z1[1, 2],subs(x1=z1[1, 1],G))));

                b := MatrixMatrixMultiply(MatrixMatrixMultiply(Transpose(evalf(subs(x2=z1[1, 2],subs(x1=z1[1, 1],G)))), evalf(subs(x2=z1[1, 2],subs(x1=z1[1, 1],Hf)))), evalf(subs(x2=z1[1, 2],subs(x1=z1[1, 1],G))));

                if b[1, 1] = 0 then

                                a[1, 1] := 0;

                else

                                a[1,1] := a[1, 1] / b[1, 1];

                end if:

                a[1, 1] := a[1, 1] / slow_constant;

                e := abs(0.0001 - evalf(subs(x2=z1[1, 2],subs(x1=z1[1, 1],obj))));

        print(e); 

                z1 := z1 - a[1, 1]*Transpose(evalf(subs(x2=z1[1, 2],subs(x1=z1[1, 1],G))));

end if:

od:

result := Vector(3);

result[1] := z1[1, 1];

result[2] := z1[1, 2];

result[3] := evalf(subs(x2=z1[1, 2],subs(x1=z1[1, 1],obj)));

ans := Minimize(obj,x1=-10..9,x2=-10..9);

[eval(f1,ans[2]), eval(f2,ans[2])];

@acer 

minmize f1 and f2 to zero is what i accept ,

but now it is a negative integer greater than 3

at least less than 1

why minimize not minimize to zero?

should i add a absolute the objective function before minimize?

@acer 

if direct search is correct and standard,

direct search 's result is what i will accept as an optimal solution

as i do not know how direct search be characterized

and i think that two objective equal to zero is what i accept

would you mind telling how direct search do?

@acer 

you are right, i am doing optimize two objective f1, f2 with 3 constraints g1,g2,g3

i haven't learnt how to determine which solutions is correct , as you can see that there are many solutions

may be this procedure should start from lambda1 = 0.1 and lambda2  = 0.1

moreover the paper do not describe when to decrease lambda1 or lambda2 during calculation,

i have an idea is that if dealing with lambda1 change to dealing with lambda2 , when dealing with lambda2 use increase first , next time chane from lambda1 to lambda2 decrease it, one increase , one decrease alternatively

however, have to know how to determine which solution is correct in order to terminate the algorithm

@acer 

is this jacobian =

[diff(diff(f, x1),x1), diff(diff(f, x1),x2), diff(diff(f, x1),x3)],

[diff(diff(f, x2),x1), diff(diff(f, x2),x2), diff(diff(f, x2),x3)],

[diff(diff(f, x3),x1), diff(diff(f, x3),x2), diff(diff(f, x3),x3)]

@Carl Love 

actually this is a question for solving eigenvectors of matrix

since sometimes after solve return 0,0,0 , i ask this question

second reason is that eigenvectors can have many kinds of solutions, some may different in sign only.

@Preben Alsholm 

actually this is a question for solving eigenvectors of matrix

since sometimes after solve return 0,0,0 , i ask this question

second reason is that eigenvectors can have many kinds of solutions, some may different in sign only.

@Markiyan Hirnyk 

steepest descent method is the standard method teaching university which is using jacobian, [df/dx1, df/x2] etc

@acer 

Hi Acer,

i got 3 different answer, one is procedure simulate by hand calculation of the paper

this procedure is not perfect, because it do not know when to decrease lambda

second is by hand trial and error

third is direct search,

how to decide which is correct?

which one is the correct answer?

with(Optimization):

f1 := -2*x1-x2; f2 := -x1-4*x2; g1 := 2*x1+3*x2-6; g2 := -x1; g3 := -x2;

penalty := lambda1*max(f1-M,0) + lambda2*max(f2-M,0) + (M^2)*(max(g1,0) + max(g2,0) + max(g3,0)):

k := 1:

s := 1:

# write into procedure

L1 := 0.5;

L2 := 0.5;

previousaction := 0; # 0 no action, 1 increase lambda1, 2 increase lambda2, 3 decrease lambda1 , 4 decrease lambda2

maxnumber := 100000;

count := 1;

succeed := 0;

delta := 0.5;

for count from 1 to maxnumber do

if succeed = 0 then

obj := eval(penalty,[lambda1=L1,lambda2=L2,M=1]);

ans := Minimize(obj,x1=-10..9,x2=-10..9);

[eval(f1,ans[2]), eval(f2,ans[2])];

if round(abs(100000*eval(f1,ans[2]))) = round(abs(100000*eval(f2,ans[2]))) then

                print("answer is");

                print(ans);

                print("L1 is ");

                print(L1);

                print("L2 is ");

                print(L2);

                print("f(x1,x2) is ");

                print(eval(f1,ans[2]), eval(f2,ans[2]));

                succeed := 1;

elif eval(f1,ans[2]) > eval(f2,ans[2]) then

                if previousaction = 2 then

                                delta := delta/2;

                                L1 := L1 - delta;

                else

                                L1 := L1 + delta;

                end if:

                previousaction := 1;

else

                if previousaction = 1 then

                                delta := delta/2;

                                L2 := L2 - delta;

                else

                                L2 := L2 + delta;

                end if:

                previousaction := 2;

end if:

end if:

od;

# by hand trial

with(VectorCalculus):

obj := eval(penalty,[lambda1=3,lambda2=0.645,M=1]);

Hf = Jacobian(Jacobian(obj, [x1, x2, x3]), [x1, x2, x3]);

Error, invalid input: VectorCalculus:-Jacobian expects its 1st argument, f,

to be of type {Vector(algebraic), list(algebraic)}, but received 3*max(0, -2*x1-x2-1)+.645*max(0, -x1-4*x2-1)+max(0, 2*x1+3*x2-6)+max(0, -x1)+max(0, -x2)

G = Matrix([[diff(obj, x1)], [diff(obj, x2)], [diff(obj, x3)]]);

ans := Minimize(obj,x1=-10..9,x2=-10..9);

[eval(f1,ans[2]), eval(f2,ans[2])];

if round(abs(100*eval(f1,ans[2]))) = round(abs(100*eval(f2,ans[2]))) then

                print("answer is");

                print(ans);

elif eval(f1,ans[2]) > eval(f2,ans[2]) then

                print("increase lambda 1");

else

                print("increase lambda 2");

end if:

# DirectSearch

with(DirectSearch):

CompromiseProgramming([f1, f2], {g1 <= 0, g2 <= 0, g3 <= 0});

@Mac Dude 

since solve can not solve (A-lambda)*x = 0 to get eigenvector

any other method to calculate eigenvector for the following matrix which contain one variable?

Matrix([[x,29.95,29.95],[29.95,29.95,0],[29.95,0,0]])

if you think it is private, please reply with email tesleft@hotmail.com

@acer 

since this post is moving backwards when more new posts,

could you reply with my email tesleft@hotmail.com?

only hope to calculate eigenvector from a matrix contain one variable, no matter the sign are correct or not