Size of SDP problem (sedumI and sostools)

[mbkv]

Hi All,

                Greetings. I have a simple question on the size of the SDP created by sostools package. (A new version has currently been released and works with the latest versions of matlab). 

The objective function used is the GoldStein-Price function (f - gamma). This is a polynomial is of 2 variables x1, x2 (n=2) and the degree of the polynomial is 2.  The demo file sosdemo3.m has been attached. I am failing to understand why the size of the SDP in this case is 226x45. 

Thanks for your help.

regards,

mbkv

[mbkv]

typo..sorry the degree of polynomial is 8.

[Johan Löfberg]

Why are you asking a sostools question on the YALMIP forum?

I'll help you by first rephrasing it is a YALMIP question

sdpvar x1 x2 gamma

f1 = x1+x2+1;

f2 = 19-14*x1+3*x1^2-14*x2+6*x1*x2+3*x2^2;

f3 = 2*x1-3*x2;

f4 = 18-32*x1+12*x1^2+48*x2-36*x1*x2+27*x2^2;

f = (1+f1^2*f2)*(30+f3^2*f4);

solvesdp(sos(f-gamma),-gamma)

What does it mean when SeDuMi says the dimension is 45 x 226. When searching for a matrix Q such that f(x)-gamma = v^T(x)Qv(x) what will the dimensions of the problem be? Since you have a degree 8 polynomial in 2 variables, v has to contain all monomials of degree 4 or less. There are  nchoosek(4+2,2) of them, i.e., 15. An 8th order polyomial in two variables has nchoosek(8+2,2)=45 terms, hence you have 45 equalities. The number 226 is 15^2+1. The 1 comes from the variable representing gamma. 

[mbkv]

I apologise for asking an sostools question here. sostools does not have a forum on its own afaik and probably the support email mentioned on their webpage is inactive.

But thank you very much for your helpful answer. 

[Johan Löfberg]

No problems, just making sure there is no misunderstanding and you think sostools is a part of YALMIP or something like that.

[mbkv]

Alright. thanks once again. 

Let me ask you some more questions to understand the problem correctly. 

So i want to transform the problem 

f(x)-gamma = v^T.Q.v

by hand into a format suitable for sedumi, namely i want to get the At, b and c. I should have an 

At of size 226x45 

b  of size 45x1

c of size  226x1

This is what should be fed to sedumi, i beleive. 

On the LHS, we have, as you explained, nchoosek(8+2, 2) gave us 45 terms. I agree with you.

What about the RHS? I symbolically multiplied 

v^T*Q*v

and i have only 29 terms. Is that how it should be done?

sedumi will have as decision variables x, where

x = [gamma q_11, q_12, .....q_15,15]

 a total of 226 decision variables

where

q_ij is the i,jth element of PSD matrix Q  

So, gamma is a free variable and rest 225 should be PSD, that is, q_ij >=0

Sedumi tries to solve: 

min c'*x

subject to A*x = b,

                   x >=0

I have two confusions now. 

1. How do i determine the entries of At(the transpose of what sedumi solves)? Should be a series of 0s,  1s and -1s when we optimize for f(x)_gamma alone?

2. How do i determine b and c?

Too newbie i know. Sorry if it was too basic. Thanks for your help.

regards,

mbkv

[Johan Löfberg]

v'*Q*v has 45 terms when seen as a polynomial in x

sdpvar x1 x2
v = monolist([x1;x2],4);
Q = sdpvar(length(v));
[coeffs,monomials] = coefficients(v'*Q*v,[x1;x2]);
sdisplay([coeffs monomials])

It is really dangerous to talk about "the number of variables" when discussing SDPs when you haven't really defined if you work in the primal or dual. When you work in the dual form, it makes sense to talk about number of variables, but in the primal, the interesting number is the number of equalities (and then you have a bunch of cones of certain dimensions which you optimize over, but you don't actually optimize over the elements in the primal cone X)

As an example. X is 1000x1000, you minimize the trace(CX) and and you only constraint is trace(HX)==1 and X>=0. A bad interpretation of this, and you might think that it is really huge since X as on the order of 1 million elements that you must optimize over. The hessian of that model would have ~10^12 elements, i.e. 1000 billions, and would not be possible to fit in anything but a supercomputer. That is the wrong to look at it though, the important thing is that you only have one constraint. The dual is max y subject to C-H*y >=0, and here we see that the problem is small also, as the number of variables only is 1. We never explicitly optimize over elements in the primal. A very very bad model would have been to see X as explicitly parameterized of half a million of variables, and extracted the data C-sum(Ai*zi) corresponding to a dual representation. Horribly intractably large.

In SOS, we typically interpret the SOS problem as a representation of a problem working in the primal. Hence, the important thing is the number of equalities connecting cone variables. In this case 45 of them, connecting a primal cone X and a small scalar cone gamma (I am assuming gamma positive for simplicity). The number 225 not that interesting. Sure, it tells us that the squared dimension of the SDP constraint 225, but complexity cubic in the size of the constraint, so reporting the square is sort of useless. It gives you the dimension to store the data, thats all.

You can tell YALMIP to derive the a variant of a dual representation instead if you which through the option sos.model and set it to 2.

solvesdp(sos(f-gamma),-gamma,sdpsettings('solver','sedumi','sos.model',2))

You will then see that you get a model with eqs m = 76, i.e, 76 equalities in the primal corresponding to 76 variables in the dual. So 76, shouldn't one get 120 variables if Q is explicitly parameterized + 1 for gamma? Yes, but then YALMIP uses the 45 equality constraint to remove that degree of freedom, so the total number of freedoms is 120+1-45 = 76. Still, the SOS model is almost twice as large when interpreted as a parameterization of the dual, instead of equalities acting on a primal cone.

[mbkv]

I greatly appreciate your detailed reply.  

1. Honestly, I am not knowledgeable about the details of interior point methods that solve SDPs. In fact, i have to consult the right  textbook to understand many terms that you have made use of in  

    in your answer. Do you have suggestions for a good textbook? 

2. Coming back again to our problem, when i talked about decision variables, i just meant that in the context of  documentation of sedumi. I am trying to get the hang of what exactly is solved by Sedumi

    and how to get data into the sedumi format. I beleive yalmip and sostools use sedumi in the backend. I am trying to understand in detail how f(x)-gamma is converted to a SDP form At,b and c so that

   it can be solved by sedumi. Essentially, this conversion is what sostools and yalmip do automatically. Why am i trying to do it manually? Because I am trying to understand an old journal paper which

   uses sedumi to perform SDP optimization and the problem is converted by hand into the sedumi format (if you are interested, please click here).

  sedumi documentation for convenience:

function [x,y,info] = sedumi(A,b,c,K,pars)

% [x,y,info] = sedumi(A,b,c,K,pars)

%

% SEDUMI  Self-Dual-Minimization/ Optimization over self-dual homogeneous

%         cones.

%

% >  X = SEDUMI(A,b,c) yields an optimal solution to the linear program

%      MINIMIZE c'*x SUCH THAT A*x = b, x >= 0

%      x is a vector of decision variables.

%      If size(A,2)==length(b), then it solves the linear program

%      MINIMIZE c'*x SUCH THAT A'*x = b, x >= 0

%

% >  [X,Y,INFO] = SEDUMI(A,b,c) also yields a vector of dual multipliers Y,

%      and a structure INFO, with the fields INFO.pinf, INFO.dinf and

%      INFO.numerr.

%

%    (1) INFO.pinf=INFO.dinf=0: x is an optimal solution (as above)

%      and y certifies optimality, viz.\ b'*y = c'*x and c - A'*y >= 0.

%      Stated otherwise, y is an optimal solution to

%      MAXIMIZE b'*y SUCH THAT c-A'*y >= 0.

%      If size(A,2)==length(b), then y solves the linear program

%      MAXIMIZE b'*y SUCH THAT c-A*y >= 0.

%

%    (2) INFO.pinf=1: there cannot be x>=0 with A*x=b, and this is certified

%      by y, viz. b'*y > 0 and A'*y <= 0. Thus y is a Farkas solution.

%

%    (3) INFO.dinf=1: there cannot be y such that c-A'*y >= 0, and this is

%      certified by x, viz. c'*x <0, A*x = 0, x >= 0. Thus x is a Farkas

%      solution.

%

%    (I)   INFO.numerr = 0: desired accuracy achieved (see PARS.eps).

%    (II)  INFO.numerr = 1: numerical problems warning. Results are accurate

%          merely to the level of PARS.bigeps.

%    (III) INFO.numerr = 2: complete failure due to numerical problems.

%

%    INFO.feasratio is the final value of the feasibility indicator. This

%    indicator converges to 1 for problems with a complementary solution, and

%    to -1 for strongly infeasible problems. If feasratio in somewhere in

%    between, the problem may be nasty (e.g. the optimum is not attained),

%    if the problem is NOT purely linear (see below). Otherwise, the reason

%    must lie in numerical problems: try to rescale the problem.

%

% >  [X,Y,INFO] = SEDUMI(A,b,0) or SEDUMI(A,b) solves the feasibility problem

%    FIND x>=0 such that A*x = b

%

% >  [X,Y,INFO] = SEDUMI(A,0,c) or SEDUMI(A,c) solves the feasibility problem

%    FIND y such that A'*y <= c

%

% >  [X,Y,INFO] = SEDUMI(A,b,c,K) instead of the constraint "x>=0", this

%      restricts x to a self-dual homogeneous cone that you describe in the

%      structure K. Up to 5 fields can be used, called K.f, K.l, K.q, K.r and

%      K.s, for Free, Linear, Quadratic, Rotated quadratic and Semi-definite.

%      In addition, there are fields K.xcomplex, K.scomplex and K.ycomplex

%      for complex-variables.

%

%    (1) K.f is the number of FREE, i.e. UNRESTRICTED primal components.

%      The dual components are restricted to be zero. E.g. if

%      K.f = 2 then x(1:2) is unrestricted, and z(1:2)=0.

%      These are ALWAYS the first components in x.

%

%    (2) K.l is the number of NONNEGATIVE components. E.g. if K.f=2, K.l=8

%      then x(3:10) >=0.

%

%    (3) K.q lists the dimensions of LORENTZ (quadratic, second-order cone)

%      constraints. E.g. if K.l=10 and K.q = [3 7] then

%          x(11) >= norm(x(12:13)),

%          x(14) >= norm(x(15:20)).

%      These components ALWAYS immediately follow the K.l nonnegative ones.

%      If the entries in A and/or c are COMPLEX, then the x-components in

%      "norm(x(#,#))" take complex-values, whenever that is beneficial.

%       Use K.ycomplex to impose constraints on the imaginary part of A*x.

%

%    (4) K.r lists the dimensions of Rotated LORENTZ

%      constraints. E.g. if K.l=10, K.q = [3 7] and K.r = [4 6], then

%          2*x(21)x(22) >= norm(x(23:24))^2,

%          2*x(25)x(26) >= norm(x(27:30))^2.

%      These components ALWAYS immediately follow the K.q ones.

%      Just as for the K.q-variables, the variables in "norm(x(#,#))" are

%      allowed to be complex, if you provide complex data. Use K.ycomplex

%      to impose constraints on the imaginary part of A*x.

%

%    (5) K.s lists the dimensions of POSITIVE SEMI-DEFINITE (PSD) constraints

%      E.g. if K.l=10, K.q = [3 7] and K.s = [4 3], then

%          mat( x(21:36),4 ) is PSD,

%          mat( x(37:45),3 ) is PSD.

%      These components are ALWAYS the last entries in x.

%

%    (a) K.xcomplex lists the components in f,l,q,r blocks that are allowed

%     to have nonzero imaginary part in the primal. For f,l blocks, these

%    (b) K.scomplex lists the PSD blocks that are Hermitian rather than

%      real symmetric.

%    (c) Use K.ycomplex to impose constraints on the imaginary part of A*x.

%

%    The dual multipliers y have analogous meaning as in the "x>=0" case,

%    except that instead of "c-A'*y>=0" resp. "-A'*y>=0", one should read that

%    c-A'*y resp. -A'*y are in the cone that is described by K.l, K.q and K.s.

%    In the above example, if z = c-A'*y and mat(z(21:36),4) is not symmetric/

%    Hermitian, then positive semi-definiteness reflects the symmetric/

%    Hermitian parts, i.e. Z + Z' is PSD.

%

%    If the model contains COMPLEX data, then you may provide a list

%    K.ycomplex, with the following meaning:

%      y(i) is complex if ismember(i,K.ycomplex)

%      y(i) is real otherwise

%    The equality constraints in the primal are then as follows:

%          A(i,:)*x = b(i)      if imag(b(i)) ~= 0 or ismember(i,K.ycomplex)

%          real(A(i,:)*x) = b(i)  otherwise.

%    Thus, equality constraints on both the real and imaginary part

%    of A(i,:)*x should be listed in the field K.ycomplex.

%

%    You may use EIGK(x,K) and EIGK(c-A'*y,K) to check that x and c-A'*y

%    are in the cone K.

%

% >  [X,Y,INFO] = SEDUMI(A,b,c,K,pars) allows you to override the default

%      parameter settings, using fields in the structure `pars'.

%

%    (1) pars.fid     By default, fid=1. If fid=0, then SeDuMi runs quietly,

%      i.e. no screen output. In general, output is written to a device or

%      file whose handle is fid. Use fopen to assign a fid to a file.

%

%    (2) pars.alg     By default, alg=2. If alg=0, then a first-order wide

%      region algorithm is used, not recommended. If alg=1, then SeDuMi uses

%      the centering-predictor-corrector algorithm with v-linearization.

%      If alg=2, then xz-linearization is used in the corrector, similar

%      to Mehrotra's algorithm. The wide-region centering-predictor-corrector

%      algorithm was proposed in Chapter 7 of

%        J.F. Sturm, Primal-Dual Interior Point Approach to Semidefinite Pro-

%        gramming, TIR 156, Thesis Publishers Amsterdam, 1997.

%

%    (3) pars.theta, pars.beta   By default, theta=0.25 and beta=0.5. These

%      are the wide region and neighborhood parameters. Valid choices are

%      0 < theta <= 1 and 0 < beta < 1. Setting theta=1 restricts the iterates

%      to follow the central path in an N_2(beta)-neighborhood.

%

%    (4) pars.stepdif, pars.w. By default, stepdif = 2 and w=[1 1]. 

%       This implements an adaptive heuristic to control ste-differentiation.

%       You can enable primal/dual step length differentiation by setting stepdif=1 or 0.

%      If so, it weights the rel. primal, dual and gap residuals as

%      w(1):w(2):1 in order to find the optimal step differentiation.

%

%    (5) pars.eps     The desired accuracy. Setting pars.eps=0 lets SeDuMi run

%      as long as it can make progress. By default eps=1e-8.

%

%    (6) pars.bigeps  In case the desired accuracy pars.eps cannot be achieved,

%     the solution is tagged as info.numerr=1 if it is accurate to pars.bigeps,

%     otherwise it yields info.numerr=2.

%

%    (7) pars.maxiter Maximum number of iterations, before termination.

%

%    (8) pars.stopat  SeDuMi enters debugging mode at the iterations specified in this vector.

%

%    (9) pars.cg      Various parameters for controling the Preconditioned conjugate

%     gradient method (CG), which is only used if results from Cholesky are inaccurate.

%    (a) cg.maxiter   Maximum number of CG-iterates (per solve). Theoretically needed

%          is |add|+2*|skip|, the number of added and skipped pivots in Cholesky.

%          (Default 49.)

%    (b) cg.restol    Terminates if residual is a "cg.restol" fraction of duality gap.

%          Should be smaller than 1 in order to make progress (default 5E-3).

%    (c) cg.refine    Number of refinement loops that are allowed. The maximum number

%          of actual CG-steps will thus be 1+(1+cg.refine)*cg.maxiter. (default 1)

%    (d) cg.stagtol  Terminates if relative function progress less than stagtol (5E-14).

%    (e) cg.qprec    Stores cg-iterates in quadruple precision if qprec=1 (default 0).

%

%    (10) pars.chol   Various parameters for controling the Cholesky solve.

%     Subfields of the structure pars.chol are:

%    (a) chol.canceltol: Rel. tolerance for detecting cancelation during Cholesky (1E-12)

%    (b) chol.maxu:   Adds to diagonal if max(abs(L(:,j))) > chol.maxu otherwise (5E5).

%    (c) chol.abstol: Skips pivots falling below abstol (1e-20).

%    (d) chol.maxuden: pivots in dense-column factorization so that these factors

%      satisfy max(abs(Lk)) <= maxuden (default 5E2).

%

%    (11) pars.vplot  If this field is 1, then SeDuMi produces a fancy

%      v-plot, for research purposes. Default: vplot = 0.

% 

%    (12) pars.errors  If this field is 1 then SeDuMi outputs some error

%    measures as defined in the Seventh DIMACS Challenge. For more details

%    see the User Guide.

%

Thanks & Regards,

mbkv

[Johan Löfberg]

1.

In this context, reading this paper is probably a good idea
http://users.isy.liu.se/en/rt/johanl/2009_OMS_DUALIZE.pdf

2.

Once you have something like

min tr C1*X1 + C2*X2 + c3*t3

tr A11*X1 + tr A12*X2 + a13*t3 = b1
tr A21*X1 + tr A22*X2 + a23*t3 = b2

where X1 and X2 are SDP cones and t3 is free variables
 
The data in sedumi is

C = [c3(:);C1(:);C2(:)]
A = [A11(:)' A12(:)' a13;A21(:)' A22(:)']
b = [b1;b2];
K.f = length(c3)
K.s = [size(A11,1) size(A12,1)]

"All" you have to figure out is what you equalities look like

trace(X1(1,1) + X2(2,3)+X2(3,2))

is [1 0 00 0 ;0 0 0 ...]*X1 + [00000;0 0 1 0 0 0;0 1 0 0 0 ]*X2 + ...

You don't want to do that for anything but trivially small problems. Why the author decide to do that in the year 2008 when there are multiple simple convenient tools for doing that escapes me.  Any small change, and they have to redo their efforts, and more advanced things like symmetry reductions, newton polytopes are missed. Sure, might want to have fast code, but then you don't use SDP+sedumi for these problem, do you? Understand small academic examples and do these by hand, and then use the tools.

[mbkv]

Your help is greatly appreciated. The pdf file you sent was very useful. I actually figured out a lot after going through the notations that clarified LP, SOCP and SDP and the nice numerical example 6.1. Now i understand that each equality constraint (obtained by equating the co-efficients of the given polynomial p(x) with the coefficients of the expanded v'*X*v) can actually be interpreted as a matrix equation  Ai * X = bi, where Ai, Xi belong to S^m and i = 1..NE, where NE is the number of equality constraints.

However, i have 2 more questions.

1. When free variables are involved, what is the relation between F^T and f which appears in the primal and dual forms?

2. sedumi returns y (as in [x,y,info] = sedumi(At, b,c)). y has been mentioned as a vector of dual variables. So how do i retreive the solution to my original problem? 

    

For instance, for the case that i study,    

The  original objective  min f(x)

                                      x

   subject to

   h(x) = 0 some equality constraints (lets say like x1+x2 = x3)

   g(x) >= 0 some inequality constraints, (for example x4>=0)

   f(x) is a polynomial function in x

Using the positivstellenstatz, we can transform that into an SDP objective:

     max       f(x) - gamma - h(x)*e(x) - g(x)*x4 is SOS

   gamma

   

   e(x) = x1 + x2

   x4 should be non-negative/g(x) is SOS

    The original problem is to solve for global x, my SDP variables are x and gamma. How do i retreive x from y returned by sedumi?

regards,

mbkv

[Johan Löfberg]

1. Skip the introduction of free variables, and see a free variable t as the difference of two positive variables, i.e., LP cones t=z1-z2, z1,z2>=0, which in the dual mean you replce Fy=f with two LP constraints. That is how it is handled internally in many solvers anyway.

2. Not possible through the sos framework (well, it is but it is messy as you have to extract a sos model using compilesos, solve the problem manually and then extract a solution from the dual of the Gramian, but then you are left with the whole theory of extrating a solution from that matrix. Not trivial). If you want to extract solutions (i.e., you are assuming the relaxation is tight and thus allows you to do this), you should use the moment relaxation framework through solvemoment instead. It implements the ideas on flat truncations etc

[mbkv]

Sorry for the late response. I had to do some groundwork before i respond to your answer. Your point 2 probably explains why the authors in that work did the extraction of the matrices A, b and c manually. The moment relaxation framework might not suit my purpose since i am looking at also having decent speed. At this point, i have actually figured out how the authors in that journal have extracted their A, b, c for input to sedumi. It was slightly tedious but i figured it out. 

I have the A, b and c for my own problem. Sedumi reports the below info:

 23 :  -1.84E+02 2.98E-13 0.000 0.1936 0.9314 0.9000   0.14  3  3  7.0E-06

 24 :  -2.09E+02 2.01E-13 0.000 0.6746 0.4672 0.4500   0.11  4  4  5.9E-06

 25 :  -3.02E+02 6.25E-14 0.000 0.3108 0.9000 0.9015  -0.09  5  5  3.8E-06

 26 :  -3.90E+02 2.93E-14 0.000 0.4689 0.9000 0.9000  -0.38 13 15  4.1E-06

 27 :  -6.13E+02 7.23E-15 0.000 0.2467 0.9055 0.9000  -0.33 15 14  2.3E-06

 28 :  -8.68E+02 2.26E-15 0.000 0.3122 0.9099 0.9000  -0.14 22 22  1.4E-06

Run into numerical problems.

iter seconds digits       c*x               b*y

 28      1.1   Inf -9.9494816005e+02 -8.6806878950e+02

|Ax-b| =   4.3e-03, [Ay-c]_+ =   8.1E-07, |x|=  5.2e+08, |y|=  2.0e+05

Detailed timing (sec)

   Pre          IPM          Post

8.154E-02    6.413E-01    5.050E-02    

Max-norms: ||b||=1.958261e+01, ||c|| = 1,

Cholesky |add|=1, |skip| = 7, ||L.L|| = 1.36086e+10.

info = 

                iter: 28

           feasratio: -0.1377

                pinf: 0

                dinf: 0

              numerr: 1

                  r0: 2.0967e-04

              timing: [0.0815 0.6413 0.0505]

             wallsec: 0.7733

              cpusec: 1.3600

    time_sedumi_call: 0.8047

So, i am now kind of confused. sedumi does not say the problem is infeasible, neither in the primal nor dual. However, i dont have my answer :) I know for sure that the objective has a global solution.solvesdp from sedumi returns me the optimal solution.

regards,

mbkv

[Johan Löfberg]

I am not sure what you mean with speed, as the moment relaxation has the same complexity as sos (unless structure is exploited in solvesos such as congruence). When I say extraction, I am not talking about extraction of the model, but extraction of a solution from a moment matrix)

SeDuMi is telling you that it cannot judge feasibility. feasratio -0.13 means "hmm, strange". 1 is fine and optimal, -1 is infeasible, in between numerical problems.

[mbkv]

"speed"? So from my understanding of moment based methods, they solve a heirarchy of convex LMI relaxations of increasing size [Globally Optimal Estimates for Geometric Reconstruction Problems]. Typically, the solution is arrived at after the second or third relaxation. So i assumed that this would mean an increase in computation time. My understanding is probably insufficient.

[Johan Löfberg]

So why would you expect sos to work then. sos and moments are two sides of the same coin. Higher order relaxations in moments means you are using higher order multipliers etc in sos

[mbkv]

1. Well, my understanding is insufficient as i told you. The sos method had worked for similar problems before. So i expected it to work for me. And i am not completely wrong.

2. I had actually messed up with the sign of the free variable (and of the equality constraints)  when assembling the At, b and c for sedumi. That is what made sedumi to report negative feasibility ratios. I corrected these problems. Now, i am able to arrive at feasratio very close to 1(1.0001 or 0.99998 etc). The solution returned is very close to the optimal values.

3. I find something strange though. The 'numerr' parameter returned by sedumi is 2 while the feasratio is 1. (The optimal solution returned was correct though) 

info = 

                iter: 25

           feasratio: 1.0000

                pinf: 0

                dinf: 0

              numerr: 2

                  r0: 0.2731

              timing: [0.0837 0.5740 0.0505]

             wallsec: 0.7082

              cpusec: 1.3300

Do you think this is normal? 

regards,

mbkv

[Johan Löfberg]

That is odd since 2 means complete failure due to numerical problems.

[mbkv]

Just to be more informative, if i run sedumi with pars.errors=1, then the code halts giving an error message about wring dimesions. This seems a bit strange to me. 

Index exceeds matrix dimensions.

Error in eigK (line 116)

    XX = x(xi+1:xi+qi);

Error in sedumi (line 797)

        info.err(2)=max(0,-min(eigK(full(x(origcoeff.K.f+1:end)),origcoeff.K)/(1+normb)));

Error in MainGlobOpti_constrained (line 29)

    [x,y,info] = sedumi(At_,b_,c_,K,pars);