You are correct, as mention in the last bullet of the Description section:

- The Minimize and Maximize commands call one of the Optimization[LPSolve], 
  Optimization[QPSolve] or Optimization[NLPSolve] commands, depending on the
  form of the input.

I should read the entire help page before complaining.

I find it interesting that method=branchandbound works here.  I would never have found that from reading the help page for Optimization[Maximize]. The Optimization,Options help page doesn't suggest that this option is available for  Minimize/Maximize, but imply it is for LPSolve and NLPSolve.

The statements in question are

y := 3*cos(4*Pi*x-1.3)+5*cos(2*Pi*x+0.5);
maximize(y, x=0..1);

I don't think a faster computer will help much.  To see this, you could do the following with command line maple (I wouldn't do this in the Standard GUI, since its output bogs down)

printlevel := 100:
maximize(y, x=0..1);

On a linux box, or with cygwin on Windows, you might try, from a shell

$ cat << EOF | maple -q | sed -n '/enter\|exit/p' | more
> y := 3*cos(4*Pi*x-1.3)+5*cos(2*Pi*x+0.5);
> printlevel := 100:
> maximize(y, x=0..1);
> EOF

Optimization:-Maximize(y,x=0..1);
                                 [2.01681169665899551, [x = 0.681550998912437866]]

While it isn't part of the comparison, we might as well use a more efficient way to generate the initial vector:

rtable(1..10\000\000, 'random'(100..1000,1));

There are two problems with (2). First, Rootof is misspelled, it should be RootOf.  Second, y is an expression, not a function (procedure); to evaluate it you should use eval.

evalf(eval(y,x=RootOf(v,0.1)));
                                 5.781103189

Not sure what you want it to do, but the type function only checks the surface type.  Use simplify:

x := sqrt(2):       
y := (x^x)^x:
simplify(y);                      2

I don't see how my approach would work with symmetries that "overlap".  It assumes that the order of application does not matter.  For your rank3 example, given, say [2,1,1], if we applied the 12 symmetry first we get [1,2,1], but then the 23 antisymmetry gives  -[1,1,2].  However, if applied in the reverse order the result is 0 (which is, I presume, what you want).

But is that a "real" symmetry?  Given say (shorter notation) we get (applying the two symmetries appropriately):

123 = 213 = -231 = -321 = 312 = 132 = -123

so any arbitrary indexes are zero.  Oh, right, that's what you meant be trivial. Regardless, I suspect this approach won't work...

Here is an alternate approach you might consider.  It assumes that you can decompose the symmetries of a tensor into group-wise symmetries and antisymmetries.  For example, for the Riemannian curvature tensor, the indices 1 and 2, as well as 3 and 4, are anti-symmetric, while the groups [12] and [34] are symmetric.  These two symmetries can be efficiently handled with Maple's builtin 'symmetric' and 'antisymmetric' index functions for tables.  Here is a rough outline of what I have in mind.

Symmetries := module()
export Anti,Symm;
local anti,symm;
    anti := table('antisymmetric'):
    symm := table('symmetric'):
    # this can be improved...
    Anti := proc( L :: seq(list), $)
    local s,S;
        S := map2(`?[]`, anti, [L]);
        if ormap(`=`, S, 0) then
            return (0,[]);
        else
            (mul(sign(s),s=S)
             , [seq(`if`(sign(s)=-1
                         , op(-s)
                         , op(s)
                        ), s=S)]
            );
        end if;
    end proc;

    Symm := proc( L :: seq(list), $)
    local s;
        [seq(op(symm[s[]]), s=[L])];
    end proc;

end module:

Now, assign a function that uses the Anti and Symm exports to handle the symmetries of the curvature tensor. This procedure returns two elements, the computed sign (-1,0,1) and a list of indices (empty if sign is 0).

Riemannian := proc(index :: list)
uses Symmetries;
local s,indx;
    (s,indx) := Anti(index[1..2],index[3..4]);
    if s = 0 then return (0,[]); end if;
    (s, map(op, Symm([indx[1..2],indx[3..4]])));
end proc:

Riemannian([1,1,2,1]);
                                     0, []
Riemannian([1,2,4,3]);
                               -1, [1, 2, 3, 4]
Riemannian([4,3,2,1]);
                                1, [1, 2, 3, 4]
Riemannian([3,4,2,1]);
                               -1, [1, 2, 3, 4]

This example was nice because the anti-symmetric indices were conveniently grouped.  However, isn't that usually the case?  One might have to do an initial permutation given an arbitrary order.

A bonus of this approach is that it will happily work with symbolic indices, however, the actual permutation generated is session dependent:

Riemannian([x,y,y,x]);
                               -1, [x, y, x, y]

Riemannian([x,y,z,z]);
                                     0, []

To check Riemannian, do

R := Array((1..4)$4, ()->Riemannian([args])): # alas, Riemannian@`[]` doesn't work
nops(map2(op,2,{entries(R)}))-1;
                                       21

which is as expected (we subtract 1 because the set included the empty list for the zero components). These symmetries don't encode the Bianchi identity, of course.

Don't pay too much attention to my last comment (which I added in an edit without too much thinking) about using sort.  That mainly works with some classes of symmetries (in particular, with completely antisymmetric tensors, but there are better ways to handle them, use the much builtin index function, antisymmetric). Note that the index table I used in the example is also sparse, so that helps somewhat.

Wouldn't there be storage issues with any approach with that sized tensor?  At least with this approach you only have to create one such large array; the actual tensor Arrays could be much smaller since all symmetries map to one location (be sure to use sparse Arrays).  You can optimize the coding.  Instead of including a scale factor (+1,-1), omit it and use a sequence to encode a positive symmetry and a list to encode an antisymmetry:

D2[1,2] := (1,2):
D2[2,1] := [1,2]:

Modify the index function accordingly. Maybe you could use a second-level of indirection to compress the index table or put some 'intelligence' into the index procedure without requiring the full computation of symmetries with each access (hand-waving).

A more efficient storage method for the indirection table is to insert +1 for positive symmetries, -1 for antsymmetries, and use the sorted indices as the base location.  That will encure a speed penalty for reading and writing, since you'll have to sort the indices.  In our example

D2[1,2] = +1
D2[2,1] = -1

I do note that your design is efficient for accessing a tensor component, since it doesn't have to use indirection as my approach does.  The piper gets paid when writing a component.  The indirect approach (my scheme) uses less storage because the Array can be sparse and all symmetries are mapped to one location.

While theoretical improvements are nice in themselves, the practical concern of whether to improve the effiicency of this routine depends on how you intend to deploy it. My assumption is that this procedure would be part of a separate procedure that is called one time, with a particular set of symmetries, and is used to create an indexing function for that class of tensor.  For example, for a tensor with four indices, you would create an index function for a 4 dimensional array with indices ranging from 1 to 4 (or similar).  To do that, you'd use this routine to partition the 4⁴ = 256 set of indices and then use that to create a special array that encodes the desired symmetry. The actual indexing function would use this array.  Because you only have to generate that array once (the indexing function can be reused for any tensor with that symmetry class), inefficiencies in it construction may not be significant. This presumes that your goal is to actually compute with the tensors, and not merely study symmetry classes of tensors.

Maybe that wasn't clear.  Here is a demo of the indexing technique I had in mind for a two-form (completely antisymmetric two-covector):

First, generate the special index table.  I do that by hand here, rather than using your code.  A zero entry (the default) indicates a zero component, otherwise the first entry in the sequence is the symmetry factor (-1,+1), the rest the indices of the stored component.

D2 := rtable(1..2,1..2,'storage=sparse'):
D2[1,2] := (1,1,2):
D2[2,1] := (-1,1,2):

 Create an index function that uses D2 (this general procedure remains the same for all tensors, you merely have to create it in a procedure that makes D2 local to it).

`index/twoform` := proc( idx :: list( posint )
                         , D :: rtable, val::list )
local d2;
    d2 := D2[idx[]];
    if _npassed = 2 then
        # retrieval
        if d2 = 0 then
            return 0;
        else
            d2[1]*D[d2[2..-1]];
        end if;
    else
        # storage
        if d2 = 0 then
            if val[] <> 0 then
                error "illegal assignment";
            else
                try
                    D[idx[]] := val[];
                catch "unable to store":
                    val[];
                end try;
            end if;
        else
            # general value
            D[d2[2..-1]] := d2[1]*val[];
        end if;
    end if;
end proc:

F := Array('twoform',1..2,1..2,'storage=sparse'):

F[1,1] := 0.;  # just checking
F[1,2] := 3:
print(F);
                                   [ 0    3]
                                   [       ]
                                   [-3    0]

See my comments on your original blog.  This is precisely what you don't want to do.  Calling permsPosNeg with each access to an index is extremely inefficient. Instead, create a special table that encodes the symmetries and do that only one time (for any symmetry class).  Use that table to create an indexing function that can be used by all tensors with that symmetry class.

While using `?[]` is neat, that it requires a second map to convert the elements to lists is an annoyance.  Consequently, I'd use seq instead.  Also, rather than creating a set of function, I'd create a function that returns a set:

(posMap,negMaps) := seq(unapply({seq([seq(x[i],i=perm)], perm in perms)},x)
                       , perms in [posPerms,negPerms]);

(posMaps,negMaps) := seq(subs('S'={seq([seq('x'[i],i=perm)], perm in perms)}, x -> S)
                         , perms in [posPerms,negPerms]);

I'm curious why you use and assign to f(t) (viz. particle(t), aysmptote(t), etc).  While it is valid Maple, it is strange in that (1) it suggests an actual function when it is not, and (2) it suggests that the quantites are functions of time, when, in fact, they are expressions equal to time (functions of x).  I might have chosen t_particle, t_asymptote, and t_photon.  Notation aside, I prefer your approach.