@brian bovril If it survives the try..catch then it augments the table of successes. If it hits an acceptable error (no solution found) then it advances to the next loop iteration. If it hits any other error then it rethrows the error so that the procedure can itself stop with that error. And it returns a sequence of the table entries.

restart:

SSP := proc(values::list)
local a, N, const, sumv, res, allsols, x;
  N := nops(values);
  const := add(x || i, i = 1 .. N);
  sumv := add(x || i*values[i], i = 1 .. N);
  for a to N do
    try
      res := Optimization:-LPSolve(sumv, {const = a, sumv = 0},
                                   assume = {binary}, depthlimit = 200);
    catch "no feasible integer point found":
      next;
    catch "no feasible point found for LP subproblem":
      next;
    catch: error;
    end try;
    allsols[a]:=res;
  end do;
  return seq([x,allsols[x]], x in indices(allsols,nolist));
end proc:

SSP([-7, -3, -2, -5, 5]);

      [2, [0, [x1 = 0, x2 = 0, x3 = 0, x4 = 1, x5 = 1]]], 

        [3, [0, [x1 = 0, x2 = 1, x3 = 1, x4 = 0, x5 = 1]]]

SSP([-7, -3, -2, 5, 8]);

       [3, [0, [x1 = 0, x2 = 1, x3 = 1, x4 = 1, x5 = 0]]]

You should contact Maple's tech support. It would probably help to tell them any language specific details, be it OS, keyboard, or Maple language pack.

@Bibinou You can suppress the system axes by not supply the axes=box option, or by supplying axes=none.

With some extra effort, you could also fake a set of axes that moved along with the rotation of the surface, by creating them as simple lines + labels in that original Q plot (and tickmarks, a little more work).

@Bibinou You can suppress the system axes by not supply the axes=box option, or by supplying axes=none.

With some extra effort, you could also fake a set of axes that moved along with the rotation of the surface, by creating them as simple lines + labels in that original Q plot (and tickmarks, a little more work).

For whatever it's worth

> f := s -> (1-s^2)/((1+s^2)*sqrt(1+s^4)):

> > value(Int(f(s), s = -1 .. 1)); # the problem (or, using just `int`)

                           1     (1/2)
                         - - Pi 2     
                           4          

> IntegrationTools:-Split(Int(f(s), s = -1 .. 1), 0);

             /  /0                          \
             | |                2           |
             | |           1 - s            |
             | |   ---------------------- ds|
             | |                    (1/2)   |
             |/-1  /     2\ /     4\        |
             \     \1 + s / \1 + s /        /

                  /  /1                          \
                  | |                2           |
                  | |           1 - s            |
                + | |   ---------------------- ds|
                  | |                    (1/2)   |
                  |/0   /     2\ /     4\        |
                  \     \1 + s / \1 + s /        /

> map(value, %);

                          1     (1/2)
                          - Pi 2     
                          4          

Whether it is practical to attempt to program such a correction is a harder queestion. Without any kind of check for validity or attempt at general correctness ( ie. just for this example!)

> q := int(f(s), s = 1 .. x);

               /        (1/2)       \                     
               |/     4\       (1/2)|                     
       1       |\1 + x /      2     |  (1/2)   1     (1/2)
     - - arctan|--------------------| 2      + - Pi 2     
       2       \        2 x         /          8
          
> d := [discont(q, x)[]];

                              [0]

d := select(z -> -1 < z and z < 1, d); # might need `is`, or might fail!

                              [0]

> value(IntegrationTools:-Split(Int(f(s), s = -1 .. 1), d));

                          1     (1/2)
                          - Pi 2     
                          4          

@Christopher2222 You asked to see that form as shown by `new` above. Something has to be done, to delay the cancellation and allow that form to be printed.

How terse, or "neat" you want to make it can depend on how much you're willing to accept specialized routines (or syntax, like with applyrule).

restart:

hold := eval(``): # or just use ``() itself (your choice)
release := z->expand(expand(numer(z)))
              /expand(expand(denom(z))):



e := (g^n) / (g^n -1);

                                n  
                               g   
                             ------
                              n    
                             g  - 1

c := g^(-n);

                              (-n)
                             g    

new := hold(numer(e)*c)/hold(denom(e)*c); 

                           / n  (-n)\   
                           \g  g    /   
                        ----------------
                        // n    \  (-n)\
                        \\g  - 1/ g    /

release(new);

                               1   
                             ------
                                 1 
                             1 - --
                                  n
                                 g 

@Christopher2222 You asked to see that form as shown by `new` above. Something has to be done, to delay the cancellation and allow that form to be printed.

How terse, or "neat" you want to make it can depend on how much you're willing to accept specialized routines (or syntax, like with applyrule).

restart:

hold := eval(``): # or just use ``() itself (your choice)
release := z->expand(expand(numer(z)))
              /expand(expand(denom(z))):



e := (g^n) / (g^n -1);

                                n  
                               g   
                             ------
                              n    
                             g  - 1

c := g^(-n);

                              (-n)
                             g    

new := hold(numer(e)*c)/hold(denom(e)*c); 

                           / n  (-n)\   
                           \g  g    /   
                        ----------------
                        // n    \  (-n)\
                        \\g  - 1/ g    /

release(new);

                               1   
                             ------
                                 1 
                             1 - --
                                  n
                                 g 

@toandhsp Yes, this demonstrates a bug in the exact symbolic `minimize` command.

You obtained a local minimum, and you did not specify a global method.

restart:	

h:=x->1+2*cos(2*x)+3*cos(4*x)+4*cos(8*x):

Optimization[Minimize](h(x),x=-Pi/12..5*Pi/12,method=branchandbound);

                 [-4.57468789382248, [x = 1.14302016898108]]

You obtained a local minimum, and you did not specify a global method.

restart:	

h:=x->1+2*cos(2*x)+3*cos(4*x)+4*cos(8*x):

Optimization[Minimize](h(x),x=-Pi/12..5*Pi/12,method=branchandbound);

                 [-4.57468789382248, [x = 1.14302016898108]]

@tringuyen Apart from using anothe name altogether, you could also use a TypeMK form of gamma.

Internally, that is a name which line-prints as `#mn(gamma)` and you can enter it as such as get the prettyprinted lowercase greek letter as output.

Another variation is `#mi("γ",fontstyle = "normal")` which prints slightly differently as output. This is what you get if you select and right-click the 2D Math input and convert it to Atomic Identifer with the context menu action (select and right-click, then) 2-D Math -> Convert To -> Atomic Identifier.

If you are daring and do not fear Maple getting confused then you can unprotect(gamma) and alter Maple's known constants and place your assumptions on `gamma`. I don't really advise doing this, however.

> restart:

> unprotect(gamma);

> constants := ({constants} minus {gamma})[];

            Catalan, FAIL, Pi, false, true, infinity

> assume(gamma<0);

> about(gamma);
Originally gamma, renamed gamma~:
  is assumed to be: RealRange(-infinity,Open(0))

>evalf(gamma);
                             gamma~

@tringuyen Apart from using anothe name altogether, you could also use a TypeMK form of gamma.

Internally, that is a name which line-prints as `#mn(gamma)` and you can enter it as such as get the prettyprinted lowercase greek letter as output.

Another variation is `#mi("γ",fontstyle = "normal")` which prints slightly differently as output. This is what you get if you select and right-click the 2D Math input and convert it to Atomic Identifer with the context menu action (select and right-click, then) 2-D Math -> Convert To -> Atomic Identifier.

If you are daring and do not fear Maple getting confused then you can unprotect(gamma) and alter Maple's known constants and place your assumptions on `gamma`. I don't really advise doing this, however.

> restart:

> unprotect(gamma);

> constants := ({constants} minus {gamma})[];

            Catalan, FAIL, Pi, false, true, infinity

> assume(gamma<0);

> about(gamma);
Originally gamma, renamed gamma~:
  is assumed to be: RealRange(-infinity,Open(0))

>evalf(gamma);
                             gamma~

@alex_01 Your (earlier, and these) plots of Y don't say anything about what the norm of XY-rhs (or Ax-y) might be.

This is not yet about your understanding of SVD. It's about your understanding of the notation and of what represents the error.

Please stop abusing symbols and notation, if possible. You originally called LeastSquares(X,Y), which means that you are trying to find a Vector T which solves X.T-Y using the least squares method. Please do not also write about mysterous Ax-y (or any other notation) as those symbols are not what your code used, and doing so only muddles everything.

Your original code has Y as the right-hand-side. Now, you have plots of estimated Y, meaning Y is what got computed (I called it T). It is very poor to first call Y the rhs of A.T=Y and then later on plot computed solution T and label it as Estimated Y.

Look, when you solve X.T=Y by least squares to obtain T then the norm of T is not by itself an indicator of how well T satisfies the equation X.T=Y. What is an somewhat ok measure is the norm of X.T-Y.

Original post: x = (A'.A)^-1.A'.y  where y is the usual rhs, A is the data matrix, and x will be the computed solution.

Code in earlier followup comment: X is data matrix, Y is the usual rhs, and data1 or data2 are the computed solutions.

Later followup Comment with plot title indicating: Y is estimated (solution).

That is three completely disparate notations you have used so far.

It is opaque what you mean by Parameter Values, in those plots. That is the first time you mentioned such. It is astounding, how often you are referring to terms which you have not overtly defined. The same goes for df (degrees of freedom -- why should everyone know that...).

Please, please, please stop using jargon from finance or statistics.

Also, please define all your terms whether from math or other disciplines, or else you risk having most of what you write remain as underfined terms awash in a sea of muddy water.

@alex_01 I'm not sure what you think that you've shown here. The normal equations can be badly conditioned enough that taking their matrix-inverse and multiplying by it can produce poor results with unnecessarily high roundoff error. The problem can also get worse, as the size increases.

I have found that some folks who resist more often the idea that it is ill-advised to numerically inverting a Matrix, just in order to then solve a linear system by multiplying with it, are those who work with computational statistics. I think it relates to some notions about covariance, and some inverse being interesing in its own right (and not for solving systems). This seems to lead to a more general belief that floating-point Matrix inversion is suitable for solving all manner of linear systems.

But the fact is that numerically inverting a Matrix and then multiplying with it, to solve a floating-point linear system, is not a good idea.

For an underdetermined linear system, problems with the normal equations can get much worse.

You should use LinearSolve instead.

I'm not sure what your plot was supposed to demonstate. You can compute forward-error, resubstituting the solution into the original X.solution -Y and taking the norm.

Do you understand that there can be some variation in a least squares solution, by minimizing the 2-norm of the resulting Vector.  The LeastSquares and the linalg:-leastsqrs commands can provide that. Did you have this in mind, perhaps, as far as what your plot might have shown!?

restart: 
randomize(): 
with(LinearAlgebra): 
with(Statistics): 

N := 10:
X := RandomMatrix(10, N, outputoptions = [datatype = float]): 
Y := Vector([seq(i, i = 1 .. 10)], datatype = float):

data1 := ( MatrixInverse(Transpose(X).X) . Transpose(X) ) . Y:
Norm(X.data1-Y);
                                         -11
                   1.30455646285554394 10   

data2 := LinearAlgebra[LeastSquares](X, Y, method=SVD):
Norm(X.data2-Y);
                                         -13
                   1.27897692436818033 10   

data3 := LinearSolve(Transpose(X).X,Transpose(X)) . Y:
Norm(X.data3-Y);
                                         -12
                   1.33582034322898835 10   


N := 20:
X := RandomMatrix(10, N, outputoptions = [datatype = float]): 
Y := Vector([seq(i, i = 1 .. 10)], datatype = float):

data1 := ( MatrixInverse(Transpose(X).X) . Transpose(X) ) . Y:
Norm(X.data1-Y);
                      75.3374023437500000

data2 := LinearAlgebra[LeastSquares](X, Y, method=SVD):
Norm(X.data2-Y);
                                         -14
                   2.30926389122032560 10   

data3 := LinearSolve(Transpose(X).X,Transpose(X)) . Y:

Norm(X.data3-Y);
                                         -15
                   8.88178419700125232 10   

Incidentally, using evalf(...,5) is also a bad idea, if all you're trying to do is round the final results to 5 decimal places. For most Maple computations, what it will do instead is cause the whole computation to use only 5 digits of working precision. That will make it much more prone to roundoff. This is not what's causing the bad result here (inverting the Matrix is). In fact you are lucky here, as UseHardware floats default value of `deduced` will still allow LinearAlgebra to use hardware double-precision. As a general rule, if you want results rounded then you should compute at usual working precision, save the result, and only as a final step apply evalf[5]. Again, this is not the major problem with what you had.

Let's move on. What transformed Matrix are you referring to? Do you understand how singular value decomposition is used in order to do least squares solving of linear systems?

@flash20001 If you have Maple 7 then your won't be able to use the legend=... option in that way. (For one thing, that typeset() facility is much more recently added functionality for Maple, and for another thing a long string in a legend will often get truncated in what is show.)

You could try textplot, and maybe pass it a string. I don't have Maple 7 to test that with, sorry. The string could be formatted using sprintf (similarly to how printf was used above) so that it contains your computed points xs and ys without your having to type in their values. I suppose that textplot & nprintf is another possible combination.