Question

I've a problem where there is a canonical solution and any rotation and mirroring (over an axis) is another solution. To avoid the problem of multiple solution due to rotation, I've set in the constraints that the vectors should be aligned with the axis when possible. The "mirror" solutions are my problem now. Basically some values can have can have a positive solution or a negative solution. This gives me 2^d solution for a problem of size d. I try to use assume and fix some values as always positive, which should solve the problem, but solve is still founding the negative solutions.

This is my code so far:

/* A parameter w0 in (0,1] */
assume(w0>=0);
assume(w0<1);
d:2$
/* d+1 extra points, the is one at x=0, for a total fo d+2 points */
s:d+1$
/* The unknown are w (weights of the first component) and x_S vectors of dimension d */
v:append([w1],makelist(x[i-s*floor(i/s)+1,floor(i/s)+1],i,0,d*s-1));
/* The different constraints */
e:append(
  /* The sum of weights is 1 */
  [w0+s*w1-1=0],
  /* Some of the components are 0 to be aligned with the axis */
  flatten(makelist(makelist(x[i,j]=0,j,i+1,d),i,1,s-1)),
  /* The mean is 0 */
  makelist(sum(x[i,j],i,1,s)=0,j,1,d),
  /* All vectors have length squared of x[1,1]^2, x[1,:] is skipped as only its first component is non-zero*/
  makelist(sum(x[i,j]^2,j,1,d)-x[1,1]^2=0,i,2,s),
  /* The diagonal of the covariance matrix is 1, w0*0 +w1*sum(x_i^2)=1*/
  makelist(w1*sum(x[i,j]^2,i,1,s)-1=0,j,1,d),
  /* The off-diagonal elements are 0. The dependancy on w1 can be eliminated since the equation is =0. */
  flatten(makelist(makelist(sum(x[i,jj]*x[i,jj+j],i,1,s)=0,j,1,d-jj),jj,1,d-1))
);
/* THIS IS NOT WORKING AS I EXPECTED, I WANT SOLUTION ONLY WITH x[i,i]>0 */
assume(x[1,1]>0,x[2,2]>0);
solution:solve(e,v)$
number_solutions=length(solution);

Is there any way to force solve to explore only some solutions of the problem?

SOLUTION: Following Robert's comment, I was able to obtain the "canonical" solution as follows:

check_canonical(sol):=block([],
  /* Extract the expression x[i,i]=... */
  diag_expr0:makelist(sublist(sol,lambda([e],(if lhs(e)=x[i,i] then true else false))),i,1,d),
  diag_expr1:flatten(diag_expr0),
  /* Get the right hand side */
  diag_expr2:makelist(rhs(diag_expr1[i]),i,1,d),
  /* Check for the rhs to be positive */
  pos_diag:sublist(diag_expr2,lambda([e],if e>0 then true else false)),
  /* If all the elelment are positive, then this is a canonical solution */
  if length(pos_diag)=d then true else false
)$
canonical_solution:flatten(sublist(solutions,check_canonical));

I'm not an expert on Maxima, but it works, although I think it would be more interesting to avoid exploring solutions that doesn't meet some criteria.

Was it helpful?

Solution

solve doesn't try to filter the set of solutions. You will have to filter the results that solve returns. Take a look at sublist.

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top