@tomleslie 

No need for DirectSearch. The min of expr^2 is of course 0 (attained for all variables 0).

@mamehrashi 

Write your problem in the text area, not in the title.
Maybe a smaller example will help.
If   P = 2*x*y + 3*x*z + 4*y*t + 5*t^2   is the polynomial in the variables x,y,z,t,
what result do you expect?

You just posted a polynomial in 16 variables, homogeneous, total degree 4,  with float coefficients.
What do you want to do with it?

@acer 

You have used the parametric version

plot3d([1,u,v],u=0..2*Pi-2*v,v=0..Pi, coords=spherical)

instead of the explicite

plot3d(1, u=0..2*Pi-2*v,v=0..Pi, coords=spherical)

It seems that you knew that the latter does not work correctly. Is it documented, or just a bug?

@Carl Love 

If we want separated colors, without interpolation, I don't think that an alternative to grid exists.

c:=piecewise(u+v<Pi, 1/6, u+v<2*Pi, 2/3, 1):
plot3d(1,u=0..2*Pi,v=0..Pi, coords=spherical, color=c, grid=[300,300], style=surface);

Compare with the default grid:

@acer Consider for simplicity  plot3d(f(x,y), x=0..1, y=0..1,  color = [c, colorscheme = [c1,c2,c3,c4]])

c = c(x,y) being the color and c1,c2,... fixed colors.

I can imagine a (useful) meaning of this for c = x (i.e. a kind of xgradient), even for c = c(x), but not for c = c(x,y).
What would be your definition? To me, it's nonsense for the moment.

plot3d( 1,u=0..2*Pi,v=0..Pi, coords=spherical, 
colorscheme= ["zgradient",[black, red, yellow, white],  colorspace="HSV"] , 
grid=[300,300]);

rr:=(x,y,z) -> piecewise(x<-1/3,0, x<1/3,1, x<=1,1):
gg:=(x,y,z) -> piecewise(x<-1/3,0, x<1/3,1, x<=1,0):
bb:=(x,y,z) -> piecewise(x<-1/3,1, x<1/3,0, x<=1,0):
plot3d(1,u=0..2*Pi,v=0..Pi, coords=spherical,
colorscheme=["xyzcoloring",[rr,gg,bb]], grid=[300,300]);

@dharr 

Yes, noticing that a "1" must exist in each row and column allows to reduce the number of the candidates.
This number can be halved (so, also the timing). And it is actually necessary to do so because otherwise each A appears twice (there are only 40320).
It is enough to change
for m to n^2 do
to

for m from k+1 to n^2 do

And preferably also k to n^2 - 1.

Edit. It's possible to reduce further the candidates by identifying two matrices A,B with A=B^+. (Maybe other symmetries too.)
BTW, for the second polynomial f:=x^8+x^4+x^3+x^2+1, A does not exist.

@Magma 

There are 2 possibilities.
1. A solution does not exist.
2. There are very few solutions and more time is needed. Actually in this case a random approach may fail, but unfortunatelly the brute force in not an option because there are > 10^11 candidates for A. [A rough estimation shows that on a regular computer almost a month could be needed].

The only solution I see is to try to reduce the number of candidates using a theoretical examination.
Or, try the "reversed" approach: take a matrix having f as characteristic polynomial (e.g. the matrix companion) and try to describe  somehow the form of its conjugates in order to find a conjugate with 10 ones.

1. Please post your code (in order to correct it). An example and the expected result would be useful too.

2. What kind of decomposition? How many terms?
What about using the elementary matrices obtained during a Gauss-Jordan elimination?

@mmcdara 

simp := eq -> `if`( eq=(0=0), undefined, eq/gcd(op(eq)));
or maybe
simp:=eq -> `if`( eq=(0=0), undefined, `if`(`*`(op(eq))=0,eq,eq/gcd(op(eq))) );

I hope that the system comes from a real world problem (otherwise just leave it).

I suggest to use DirectSearch (SolveEquations). Include the domain where you are expecting a solution (something like mu1 = 0.1 .. 20, ...), set Digits (>=15)  and let it go. Probably in a few days a solution will appear.

@acer Thanks for the answer. The second question should have not been posed because obviously Seq only has to deal with unevaluted functions and plot is not involved.
 

What could be an explanation for the following two facts?

1. On my computer, I get:

Error, (in plot) attempting to assign to `ArrayTools:-Alias` which is protected

Executing frames := ...  again, the error disappears. The artifacts are then present.

2. If I replace 
frames := Threads:-Seq('frame'(t), t=0..1, 0.05):

everything works fine, no artifacts. (I don't know whether Seq has special evaluation rules. Maybe in this case it is seq which is actually executed.).

@minhthien2016 
If A is a column vector, A^+  or  A^%T  is the row vector  (i.e. "transposed").
So, 

(P-H)^+ . (B-A)

is the (real) dot product (bilinear form).