Coding

ClearAll["Global`*"]

\[Phi]0[x_] := 1/Sqrt[2 \[Pi]]*E^(-(x^2/2));

\[Phi][x_] :=

1/(Sqrt[2 \[Pi]]*\[Sigma])*E^(-((x - \[Mu])^2/(2*\[Sigma]^2)));

HermiteM[x_, n_] := Expand[2^(-n/2)*HermiteH[n, x/Sqrt[2]]];

GCPolynomial[x_, n_] :=

1 + Sum[Subscript[\[Gamma], k]/k!*HermiteM[x, k], {k, 3, n}];

order = 16;

GCPolynomial[x, order];

Coefflist = CoefficientList[GCPolynomial[x, order], x];

ParaGamma = Table[Subscript[\[Gamma], k], {k, 3, order}];

n = 0;

matrice = {};

While[n <= order/2,

matrice =

Append[matrice,

Join[Table[Subscript[q, m, n], {m, 0, n - 1}],

Table[Subscript[q, n, m], {m, n, order/2}]]];

n++

];

Clear[n];

(*matrice//MatrixForm*)

Flamatrice = Flatten[matrice];

matriceCoeff =

Table[Sum[

Select[Flamatrice, #[[2]] + #[[3]] == n &][[m]], {m, 1,

Length[Select[Flamatrice, #[[2]] + #[[3]] == n &]]}], {n, 0,

order}];

(*matriceCoeff//TableForm*)

relation = {};

nn = 1;

relation = Append[relation, matriceCoeff[[nn]] -> Coefflist[[nn]]];

relation =

Append[relation,

matriceCoeff[[Length[Coefflist] - nn + 1]] ->

Coefflist[[Length[Coefflist] - nn + 1]]];

nn = 2;

relation =

Append[relation, 1/2*matriceCoeff[[nn]] -> 1/2*Coefflist[[nn]]];

relation =

Append[relation,

1/2*matriceCoeff[[Length[Coefflist] - nn + 1]] ->

1/2*Coefflist[[Length[Coefflist] - nn + 1]]];

ParaLambda = {};(*用来收集一共用了多少Lambda*)

For[n = 3, n <= 2 + IntegerPart[((Length[Coefflist] - 1)/2 - 2)/2],

n++,

ResParaTotal =

Table[

Subscript[\[Lambda], m], {m, 1 + 2*(n - 3)*(n - 2),

2*(n - 3)*(n - 2) + 4*(n - 2)}];

ParaLambda = Join[ParaLambda, ResParaTotal];

ResPara = Take[ResParaTotal, Length[ResParaTotal]/2];

ResPara1 = Take[ResParaTotal, -Length[ResParaTotal]/2];

(*2 Subscript[q, 0,4]+2 Subscript[q, 1,3]+Subscript[q, 2,2],

最后一项系数不是2的情形*)

nn = 2*(n - 1) - 1;

length = Length[matriceCoeff[[nn]]];

(*这一条对角线上有多少个q等着去分配*)

relation =

Append[relation, matriceCoeff[[nn]][[length]] -> 2 ResPara[[1]]];

(*先把最后一个，没有系数2的那个q，即对角线上的元素分配一个\[Lambda]，选择ResPara中的第一个\[Lambda]*)

If[length >= 3(*至少要有三项，才至少有一个给中间的配\[Lambda]*),

relation =

Join[relation,

Table[

1/2*matriceCoeff[[nn]][[m]] -> ResPara[[m]], {m, 2,

Length[matriceCoeff[[nn]]] - 1}]];

,];

(*除掉第一个和最后一个q，如果中间还有元素则逐个分配\[Lambda]，选择ResPara中的第二个\[Lambda]开始一直到*)

\

relation =

Append[relation,

1/2*matriceCoeff[[nn]][[1]] ->

1/2*Coefflist[[nn]] - Sum[ResPara[[m]], {m, 1, length - 1}]];

(*把第一个q分配给GC的系数，再减去前面所有的\[Lambda]*)

(*对应的倒数第三也需要处理，因为是对称的处理。*)

relation =

Append[relation,

matriceCoeff[[Length[Coefflist] - nn + 1]][[length]] ->

2 ResPara[[Length[matriceCoeff[[nn]]]]]];

(*先把最后一个，没有系数2的那个q，即对角线上的元素分配一个\[Lambda]，选择ResPara中的第一个\[Lambda]*)

If[length >= 3(*至少要有三项，才至少有一个给中间的配\[Lambda]*),

relation =

Join[relation,

Table[

1/2*matriceCoeff[[Length[Coefflist] - nn + 1]][[m]] ->

ResPara[[m + Length[matriceCoeff[[nn]]] - 1]], {m, 2,

Length[matriceCoeff[[nn]]] - 1}]];

,];

(*除掉第一个和最后一个q，如果中间还有元素则逐个分配\[Lambda]，选择ResPara中的第二个\[Lambda]开始一直到*)

\

relation =

Append[relation,

1/2*matriceCoeff[[Length[Coefflist] - nn + 1]][[1]] ->

1/2*Coefflist[[Length[Coefflist] - nn + 1]] -

Sum[ResPara[[m]], {m, length, Length[ResPara]}]];

(*把第一个q分配给GC的系数，再减去前面所有的\[Lambda]*)

(*2 Subscript[q, 0,5]+2 Subscript[q, 1,4]+2 Subscript[q, 2,3],

系数都是2的情形*)

nn = 2*(n - 1);

length = Length[matriceCoeff[[nn]]];

(*这一条对角线上有多少个q等着去分配*)

relation =

Append[relation, 1/2*matriceCoeff[[nn]][[length]] -> ResPara1[[1]]];

(*先把最后一个，没有系数2的那个q，即对角线上的元素分配一个\[Lambda]，选择ResPara中的第一个\[Lambda]*)

If[length >= 3(*至少要有三项，才至少有一个给中间的配\[Lambda]*),

relation =

Join[relation,

Table[

1/2*matriceCoeff[[nn]][[m]] -> ResPara1[[m]], {m, 2,

Length[matriceCoeff[[nn]]] - 1}]];

,];

(*除掉第一个和最后一个q，如果中间还有元素则逐个分配\[Lambda]，选择ResPara中的第二个\[Lambda]开始一直到*)

\

relation =

Append[relation,

1/2*matriceCoeff[[nn]][[1]] ->

1/2*Coefflist[[nn]] - Sum[ResPara1[[m]], {m, 1, length - 1}]];

(*把第一个q分配给GC的系数，再减去前面所有的\[Lambda]*)

(*对应的倒数第三也需要处理，因为是对称的处理。*)

relation =

Append[relation,

1/2*matriceCoeff[[Length[Coefflist] - nn + 1]][[length]] ->

ResPara1[[Length[matriceCoeff[[nn]]]]]];

(*先把最后一个，没有系数2的那个q，即对角线上的元素分配一个\[Lambda]，选择ResPara中的第一个\[Lambda]*)

If[length >= 3(*至少要有三项，才至少有一个给中间的配\[Lambda]*),

relation =

Join[relation,

Table[

1/2*matriceCoeff[[Length[Coefflist] - nn + 1]][[m]] ->

ResPara1[[m + Length[matriceCoeff[[nn]]] - 1]], {m, 2,

Length[matriceCoeff[[nn]]] - 1}]];

,];

(*除掉第一个和最后一个q，如果中间还有元素则逐个分配\[Lambda]，选择ResPara中的第二个\[Lambda]开始一直到*)

\

relation =

Append[relation,

1/2*matriceCoeff[[Length[Coefflist] - nn + 1]][[1]] ->

1/2*Coefflist[[Length[Coefflist] - nn + 1]] -

Sum[ResPara1[[m]], {m, length, Length[ResPara1]}]];

(*把第一个q分配给GC的系数，再减去前面所有的\[Lambda]*)

];

(*下面处理最后一个对角线上元素*)

If[OddQ[order/2],

ResParaTotal =

Table[

Subscript[\[Lambda], m], {m, 1 + 2*(n - 3)*(n - 2),

2*(n - 3)*(n - 2) + 4*(n - 2)}];

ResPara = Take[ResParaTotal, Length[ResParaTotal]/2];

ResPara1 =

Take[ResParaTotal, {Length[ResParaTotal]/2 + 1, (

3*Length[ResParaTotal])/4}];

ParaLambda = Join[ParaLambda, ResPara];

ParaLambda = Join[ParaLambda, ResPara1];

nn = 2*(n - 1) - 1;

length = Length[matriceCoeff[[nn]]];

(*这一条对角线上有多少个q等着去分配*)

relation =

Append[relation, matriceCoeff[[nn]][[length]] -> 2 ResPara[[1]]];

(*先把最后一个，没有系数2的那个q，即对角线上的元素分配一个\[Lambda]，选择ResPara中的第一个\[Lambda]*)

If[length >= 3(*至少要有三项，才至少有一个给中间的配\[Lambda]*),

relation =

Join[relation,

Table[

1/2*matriceCoeff[[nn]][[m]] -> ResPara[[m]], {m, 2,

Length[matriceCoeff[[nn]]] - 1}]];

,];

(*除掉第一个和最后一个q，如果中间还有元素则逐个分配\[Lambda]，选择ResPara中的第二个\[Lambda]开始一直到*)

\

relation =

Append[relation,

1/2*matriceCoeff[[nn]][[1]] ->

1/2*Coefflist[[nn]] - Sum[ResPara[[m]], {m, 1, length - 1}]];

(*把第一个q分配给GC的系数，再减去前面所有的\[Lambda]*)

(*对应的倒数第三也需要处理，因为是对称的处理。*)

relation =

Append[relation,

matriceCoeff[[Length[Coefflist] - nn + 1]][[length]] ->

2 ResPara[[Length[matriceCoeff[[nn]]]]]];

(*先把最后一个，没有系数2的那个q，即对角线上的元素分配一个\[Lambda]，选择ResPara中的第一个\[Lambda]*)

If[length >= 3(*至少要有三项，才至少有一个给中间的配\[Lambda]*),

relation =

Join[relation,

Table[

1/2*matriceCoeff[[Length[Coefflist] - nn + 1]][[m]] ->

ResPara[[m + Length[matriceCoeff[[nn]]] - 1]], {m, 2,

Length[matriceCoeff[[nn]]] - 1}]];

,];

(*除掉第一个和最后一个q，如果中间还有元素则逐个分配\[Lambda]，选择ResPara中的第二个\[Lambda]开始一直到*)

\

relation =

Append[relation,

1/2*matriceCoeff[[Length[Coefflist] - nn + 1]][[1]] ->

1/2*Coefflist[[Length[Coefflist] - nn + 1]] -

Sum[ResPara[[m]], {m, length, Length[ResPara]}]];

(*把第一个q分配给GC的系数，再减去前面所有的\[Lambda]*)

(*下面处理对角线上元素！！！！！！！！！！！！！！！！！*)

nn = 2*(n - 1);

length = Length[matriceCoeff[[nn]]];

(*这一条对角线上有多少个q等着去分配*)

relation =

Append[relation, 1/2*matriceCoeff[[nn]][[length]] -> ResPara1[[1]]];

(*先把最后一个，没有系数2的那个q，即对角线上的元素分配一个\[Lambda]，选择ResPara中的第一个\[Lambda]*)

If[length >= 3(*至少要有三项，才至少有一个给中间的配\[Lambda]*),

relation =

Join[relation,

Table[

1/2*matriceCoeff[[nn]][[m]] -> ResPara1[[m]], {m, 2,

Length[matriceCoeff[[nn]]] - 1}]];

,];

(*除掉第一个和最后一个q，如果中间还有元素则逐个分配\[Lambda]，选择ResPara中的第二个\[Lambda]开始一直到*)

\

relation =

Append[relation,

1/2*matriceCoeff[[nn]][[1]] ->

1/2*Coefflist[[nn]] - Sum[ResPara1[[m]], {m, 1, length - 1}]];

(*把第一个q分配给GC的系数，再减去前面所有的\[Lambda]*)

,(*下面是OddQ[order/2]=False的情况*)

(*下面"只"处理对角线上元素！！！！！！！！！！！！！！！！！*)

nn = 2*(n - 1) - 1;

ResParaTotal =

Table[

Subscript[\[Lambda], m], {m, 1 + 2*(n - 3)*(n - 2),

2*(n - 3)*(n - 2) + 4*(n - 2)}];

ResPara = Take[ResParaTotal, Length[ResParaTotal]/4];

ParaLambda = Join[ParaLambda, ResPara];

nn = 2*(n - 1) - 1;

length = Length[matriceCoeff[[nn]]];

(*这一条对角线上有多少个q等着去分配*)

relation =

Append[relation, matriceCoeff[[nn]][[length]] -> 2 ResPara[[1]]];

(*先把最后一个，没有系数2的那个q，即对角线上的元素分配一个\[Lambda]，选择ResPara中的第一个\[Lambda]*)

If[length >= 3(*至少要有三项，才至少有一个给中间的配\[Lambda]*),

relation =

Join[relation,

Table[

1/2*matriceCoeff[[nn]][[m]] -> ResPara[[m]], {m, 2,

Length[matriceCoeff[[nn]]] - 1}]];

,];

(*除掉第一个和最后一个q，如果中间还有元素则逐个分配\[Lambda]，选择ResPara中的第二个\[Lambda]开始一直到*)

\

relation =

Append[relation,

1/2*matriceCoeff[[nn]][[1]] ->

1/2*Coefflist[[nn]] - Sum[ResPara[[m]], {m, 1, length - 1}]];

(*把第一个q分配给GC的系数，再减去前面所有的\[Lambda]*)

];

Clear[n];

matriceNew = matrice /. relation // Expand;

matriceNew // MatrixForm

(*做检验，生成的矩阵，再让其对角线相加，看看是不是等于GC的各阶系数*)

n = 1;

list = {};

While[n <= Length[Coefflist],

If[n <= (Length[Coefflist] + 1)/2,

list = Append[list, Sum[matriceNew[[n - m + 1]][[m]], {m, 1, n}]];

,

list =

Append[list,

Sum[

matriceNew[[m]][[n - m + 1]], {m,

n - (Length[Coefflist] + 1)/2 + 1, (Length[Coefflist] + 1)/

2}]];

];

n++];

Clear[n];

(*list==Coefflist(*判断是不是等于原GC的各阶系数*)*)

(*SymmetricMatrixQ[matriceNew](*判断是不是对称矩阵*)*)

Print["All Parameters =", Join[ParaGamma, ParaLambda]];

Print["分解矩阵为："];

Subscript[matrice, 0] =

matriceNew /.

Join[Table[

Subscript[\[Lambda], n] -> 0, {n, 1, Length[ParaLambda]}],

Table[Subscript[\[Gamma], n] -> 0, {n, 3, Length[ParaGamma] + 2}]];

(*Subscript[matrice,0]//MatrixForm*)

k = 3;

While[k <= Length[ParaGamma] + 2,

Subscript[matriceGammaPre, k] =

matriceNew - Subscript[matrice, 0] /.

Join[

Table[Subscript[\[Lambda], n] -> 0, {n, 1, Length[ParaLambda]}],

Join[Table[Subscript[\[Gamma], n] -> 0, {n, 3, k - 1}],

Table[

Subscript[\[Gamma], n] -> 0, {n, k + 1,

Length[ParaGamma] + 2}]]];

Subscript[matriceGamma, k] =

Subscript[matriceGammaPre, k]/Subscript[\[Gamma], k];

(*Print[Subscript[matriceGamma,k]//MatrixForm];*)

k++];

k = 1;

While[k <= Length[ParaLambda],

Subscript[matriceLambdaPre, k] =

matriceNew - Subscript[matrice, 0] /.

Join[

Table[

Subscript[\[Gamma], n] -> 0, {n, 3, Length[ParaGamma] + 2}],

Join[Table[Subscript[\[Lambda], n] -> 0, {n, 1, k - 1}],

Table[

Subscript[\[Lambda], n] -> 0, {n, k + 1, Length[ParaLambda]}]]];

Subscript[matriceLambda, k] =

Subscript[matriceLambdaPre, k]/Subscript[\[Lambda], k];

(*Print[Subscript[matriceLambda,k]//MatrixForm];*)

k++];

Clear[k];


res = SemidefiniteOptimization[Subscript[\[Gamma], 3],

VectorGreaterEqual[{Activate[Subscript[matrice, 0] \!\(\*

TagBox["+",

"InactiveToken",

BaseStyle->"Inactive",

SyntaxForm->"+"]\)

Join[ParaLambda, ParaGamma] .

Join[

Table[

Subscript[matriceLambda, k], {k, 1, Length[ParaLambda]}],

Table[

Subscript[matriceGamma, k], {k, 3, Length[ParaGamma] + 2}]],

Unevaluated[Plus]], 0}, {"SemidefiniteCone", (

Length[Coefflist] + 1)/2}], Join[ParaLambda, ParaGamma],

Method -> {"CSDP", Tolerance -> 0.041, MaxIterations -> 1000000}]

Print["最小的 Skewness = "];

Subscript[\[Gamma], 3] /. res