@vv: I guess you have a point. I did not notice that at all when reading the question.

@mskalsi: Sorry, for first returning back now. I am afraid that I have no idea of how to proceed on this problem.

@Preben Alsholm: Thanks for your great effort. You have given me some new avenues to pursue, not the least the function evalindets, which I had all forgotten about. I will shut my computer down for today, but return tomorrow to look further into the matters, and hopefully report back on some progress.

PS: You suspect correctly: I have several functions. In fact I have eight functions and four operators, the latter of which can be combined in 16 ways, thus yielding a total of 8 x 16 = 128 substitutional relations to be taken care of.

Update (reporting back): I have never used element-wise (~) operators before, so it took me a little time to fully decipher your latter example. However, your code works, but as far as I can tell only for a single operator at a time. This can readily be remedied, though, the following being an example thereof:

opers := {tau,chi[1],chi[2],chi[3]};   # An example with four operators
funcs := {psi1,psi2,psi3,psi4};        # An example with four functions
conv := proc(expr)
   local types,Subs,g,res;
   types := map(x -> identical(x),opers union funcs);
   Subs := `union`(
      {seq(opers[1]*func = func(-t,x,y,z),func in funcs)},   # t -> -t
      {seq(opers[2]*func = func(t,-x,y,z),func in funcs)},   # x -> -x
      {seq(opers[3]*func = func(t,x,-y,z),func in funcs)},   # y -> -y
      {seq(opers[4]*func = func(t,x,y,-z),func in funcs)}    # z -> -z
   );
   res := evalindets(expr,`*`,s -> g(selectremove(type,s,types)));
   res := eval(res,Subs);
   eval(%,g = `*`);
end proc:
expr := [I*f(x)*tau*psi1,h(t)*chi[2]*psi2];
exprConv := conv(expr);

I have taken the liberty to rewrite some of your code in a way that I am more accustomed to. The last three lines of your code, which are the really clever and critical ones, I have left essentially untouched. Introducing the auxiliary g, and latter replace it with `*`, would never have occured to me; in fact, I would not even have considered it possible. Anyway, thanks a lot for your help. It is a shame that your solution figures in a reply to which I am not able to give a thumbs-up point (which it certainly deserves).

A little note: The elements of funcs above are given without any arguments; they are given in the abstract representation, I use. By mistake, the original question of mine was posed differently, with functions with arguments. However, it makes no difference to the issue or its solution.

@vv: I am not at all qualified to judge on these matters of logical versus internal representations. And that is not, of course, meant reproachfully.

@Preben Alsholm: The example I show is only a toy example meant to flesh out the issue. But * is multiplication also in my actual case. What I am doing, is setting up a representation for operators acting on some functions. Besides the time reversal operator, I also have space inversion operators.

For long stretches of calculation, it is far easier to work with such a representation, than to work with real operators: in the representation, composition of operators is just multiplication of their symbolic placeholders (because the operators all commute). At the end, I would then like to convert back to non-symbolic, real expressions.

@vv: Interesting!

"The unpleasant fact is that op and nops do not reflect in this case the internal representation": Is that as it should be, or could it be considered a bug?

@Preben Alsholm: So it seems that I am trying to do something with eval (or subs) that it is not really intended to be able to do. It would then seem that algsubs is the only command to use (or am I mistaken?), but as I have written in an additional note to the question posed, I would like to avoid that because it performs only a single substitution at a time, and I have more than a hundred of them; having to loop through them would seem to me quite extraordinary.

@Preben Alsholm: Thanks for your reply.

Using 3 instead of I is covered by "Note that if I is replaced by any real number, then the substitution works quite as expected" in my question. So the fact that that works does not surprise me.

Your suggestion dividing by tau does not really work, as the substitution should have no effect on f(t) in isolation.

Could you explain what 'optimal' refers to?

@vv: I think you have nailed it. The solution clearly gives seven of the eight functions in terms of operators acting on the remaining one, in the specific case above, f1,f2,f3,f4,f5,f6,f7 in terms of f8. I most certainly have used Groebner wrong. I have had zero experience with Groebner until today.

Something remains confusing to me, though. When trying to translate your choice of variable names into the ones I use, I do not get a corresponding result: the method does not solve for seven of the functions in terms of the remaining one. Quite plausibly, though, it is just me fooling around. I will return, if I cannot find the bug in my thinking.

Update I: I have located the 'bug' now. I was not aware that plex() was associated with lexicographic ordering; it first became clear to me when reading the help page. If just putting indets(EQS) inside plex(), I get the wrong ordering, with operators coming before the functions, rather than functions coming before operators, as it has to be (although at present I have no clear understanding of why this has to be so).

Update II: It suddenly dawned on me that if I wanted to express, say, f2,f3,f4,f5,f6,f7,f8 in terms of f1, rather than f1,f2,f3,f4,f5,f6,f7 in terms of f8, then probably I just had to put f1 last in the lexicographic ordering of the eight functions. And it works, giving me even more control. Nice! Reassuringly, it yields the manually obtained solution given elsewhere in this thread.

@Markiyan Hirnyk: Thanks for keep on trying. Unfortunately, neither the full parametric method, as you link to, or the eliminate method above, avoids assigning to the operators (the Omega's) something: For the former method, it assigns +1 or -1; for the latter method, it even assigns fractions of some of the functions. Neither assignment makes sense. The operators should be left entirely unassigned (as they are in the manually obtained explicit solution I gave above).

@Markiyan Hirnyk: No, why should it be empty talk? The following code (which can be appended to the above given worksheet) verifies my claim:

sol := {
   f[2] = Omega[3]*Omega[4]*f[1],
   f[3] = Omega[2]*Omega[3]*Omega[4]*f[1],
   f[4] = Omega[2]*f[1],
   f[5] = Omega[1]*Omega[2]*Omega[4]*f[1],
   f[6] = Omega[1]*Omega[2]*Omega[3]*f[1],
   f[7] = Omega[1]*Omega[3]*f[1],
   f[8] = Omega[1]*Omega[4]*f[1]
};
sideRels := {seq(Omega[i]^2 = 1,i = 1..4)};
simplify(eval(EQS,sol),sideRels);

@Markiyan Hirnyk: Here they are:

SolveWithSideRelations.mw

SolveWithSideRelations.txt

Although not being my original code, the worksheet can stand alone, and it produces the system of equations at hand.

@Markiyan Hirnyk: Ok, here they come:

The Omega's are the four mutually commuting operators, each with square equal to the identity operator; and the f's are the eight functions on which they act. As earlier expressed, what I want is to have the f's given in terms of one another, as say f7 = Omega1*Omega3*f1, etc. Actually, having solved the equations manually, I know that, say, f2,f3,f4,f5,f6,f7,f8 can all be expressed in terms of various combinations of the four operators acting on f1. So the solution has only a single functional degree of freedom.

@vv: Thanks for your second try with the Groebner bases. Unfortunately, when appropiately edited to the actual situation I have, your code returns only the third argument list of indeterminates, i.e., in my case the names of the eight functions (on which the four operators are supposed to act).