@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).

@tomleslie 

The big problem is that the OP has copied the 2D math as text and has pasted here.
(In my opinion the 2D math is more difficult to be used correctly for a beginner than the old 1D; it should not be the default!).

BTW, in 2D,  f(x) := ...   is converted into an arrow procedure definition f := x -> ...  Of course f(x+1) := ... is superfluous and represents an useless entry in the remember table.

@dharr 

Yes, it's strange. It seems that when the point is infinity, taylor calls `series/infinity` and does not check the polynomiality of the result.
IIRC, in the past taylor was a regular procedure but now it is builtin and it's not posible to see its code.
Pi versus pi is not relevant here:  taylor(x^(1/x), x=infinity);

Actually we can obtain a little more general result:

If P is a polynomial in x with deg(P,x) = 2*n  then it has a unique representation
         P = a*Q^2 + R, 
where a is constant, Q, R are polynomials, Q is monic, deg(Q,x) =n, deg(R,x) <= n-1.

The procedure CompleteSquare rejects the representation if R is not constant.

SQR1:=proc(P::polynom(anything,x), x::name)
local n:=degree(P,x)/2, q,r,Q,R,a;
if not(type(n,posint)) then error "degree(P) must be even" fi;
a:=lcoeff(P,x);
Q:=add(q[k]*x^k,k=0..n-1) + x^n;
R:=add(r[k]*x^k,k=0..n-1);  #q,r ==> 2*n unknowns
solve({coeffs(expand(P-a*Q^2-R),x)}, {seq(q[k],k=0..n-1),seq(r[k],k=0..n-1)});
a*eval(Q,%)^2+eval(R,%)
end:

SQR1(2*x^6-8*x^4-5*x^3-9*x^2-4*x-8,x);
SQR1( a*x^2+b*x+c, x );
SQR1(54*x^6-216*x^5-72*x^4+684*x^3+168*x^2-288*x+51, x);
               

@Vitreg2 

The taylor command does not (and cannot) return such things.
Please post your worksheet (use the green arrow on the toolbar).

@Kitonum 

Yes, but e.g. if you want to plot3d it, you must use first expand  etc.

@Magma 

Inside the iteration you can do anything with v which contains the list [i1,i2,...]. In CartesianProduct the only difference is that v is an Array.

f := (i, r) -> add(i[j]*r[j], j=1..nops(r));

Q:=proc(k::list(nonnegint), r::list)
local P,i,v;
P := combinat:-cartprod([seq([seq(0..i)],i=k)]);
while not P[finished] do 
  v :=  P[nextvalue]();
  print( f(v,r) ); 
end do;
end:

Q([2,3],[r1,r2]);
                               0
                               r2
                              2 r2
                              3 r2
                               r1
                            r1 + r2
                           r1 + 2 r2
                           r1 + 3 r2
                              2 r1
                           2 r1 + r2
                          2 r1 + 2 r2
                          2 r1 + 3 r2

@Joe Riel 
This works:

r2_r1 := r -> rhs(Optimization:-Maximize(eval(Ept,r1=r), r2=1000..1700)[2][]):
plot(r2_r1, 1000..1700);