[R] Quadratic programming

[Maija Sirkjärvi]

Thank you for giving me your time!

The problem is the quadratic optimization part. Something goes wrong along
the way. In C++ loops run from 0 and in R they run from 1, and I've tried
to take that into account. Still I'm having hard time figuring out the
mistake I make, cause I get a result from my R code. It's just not the same
that I get with the C++.

Here are the quadratic optimization parts for both codes.

C++

  /* Set Up Quadratic Programing Problem */
Vector<double> hSmooth(J);
for(int j=0; j<J; j++) hSmooth(j) = -pow(kr.Phi(j),Rho);
Vector<double> Q(J);
for(int j=0; j<J; j++) Q(j) = exp(-Rho * (Beta * pow(Price(j),Eta + One) -
One) / (One + Eta));
SymmetricMatrix<double> H(J,Zero);
Vector<double> c(J,Zero);
Matrix<double> Aeq(0,J);
Vector<double> beq(0);
Matrix<double> Aneq(2*J-3,J,Zero);
Vector<double> bneq(2*J-3);
Vector<double> lb(J,-Inf);
Vector<double> ub(J,Inf);
for(int j=0; j<J; j++) H(j,j) = One;
for(int j=0; j<J; j++) c(j) = -hSmooth(j);

for(int j=1; j<J; j++)
{
Aneq(j-1,j-1) = -One;
Aneq(j-1,j)   = One;
bneq[j-1]     = Delta1;
}
for(int j=2; j<J; j++)
{
Aneq(J-1+j-2,j)   = -One / (Q(j) - Q(j-1));
Aneq(J-1+j-2,j-1) = One / (Q(j) - Q(j-1)) + One / (Q(j-1) - Q(j-2));
Aneq(J-1+j-2,j-2) = -One / (Q(j-1) - Q(j-2));
bneq[J-1+j-2]     = Delta2;
}

/* Solve Constrained Optimization Problem Using Quadratic Programming */
MinQuadProg qp(c,H,Aeq,beq,Aneq,bneq,lb,ub);
qp.PrintLevel = 0;
qp.Solve();

And R

J <- length(Price)
hs <- numeric(J)
for(j in 1:J){
  hs[j] <-(-(gEst$KernelRegPartLin..Phi[j]^(-0.1)))
}
hs

Q <- rep(0,J)
for(j in 1:(length(Price))){
  Q[j] <- exp((-0.1) * (Beta *Price[j]^(Eta + 1) - 1) / (1 + Eta))
}
Q
plot(Q)
Dmat <- matrix(0,nrow= J, ncol=J)
diag(Dmat) <- 1
dvec <- -hs
Aeq <- 0
beq <- 0
Amat <- matrix(0,J,2*J-3)
bvec <- rep(0,2*J-3)

for(j in 2:nrow(Amat)){
  Amat[j-1,j-1] = -1
  Amat[j,j-1] = 1
}

for(j in 3:nrow(Amat)){
  Amat[j,J+j-3] = -1/(Q[j]-Q[j-1])
  Amat[j-1,J+j-3] = 1/(Q[j]-Q[j-1])
  Amat[j-2,J+j-3] = -1/(Q[j-1]-Q[j-2])
}

for(j in 2:nrow(bvec)) {
  bvec[j-1] = Delta1
}
for(j in 3:nrow(bvec)) {
  bvec[J-1+j-2] = Delta2
}

solution <- solve.QP(Dmat,dvec,Amat,bvec)





ke 23. syysk. 2020 klo 9.12 Abby Spurdle (spurdle.a using gmail.com) kirjoitti:

> > I'm trying to replicate a C++ code with R.
>
> Notes:
> (1) I'd recommend you make the code more modular.
> i.e. One function for initial data prep/modelling, one function for
> setting up and solving the QP, etc.
> This should be easier to debug.
> (However, you would probably have to do it to the C++ code first).
> (2) Your R code is not completely reproducible.
> i.e. AEJData
> (3) For the purposes of a reproducible example, your code can be
> simplified.
> i.e. Only one contributed R package should be attached.
>
> Regardless of (1) above, you should be able to identify at what point
> the C++ and R code becomes inconsistent.
> The simplest approach is to add print-based functions into both the
> C++ and R code, and print out state data, at each major step.
> Then all you need to do is compare the output for both.
>
> > Is there a better/another package for these types of problems?
>
> I'm not sure.
> After a quick search, this is the best I found:
>
> scam::scam
> scam::shape.constrained.smooth.terms
>

	[[alternative HTML version deleted]]