Thank you both for your responses!
For practical purposes I thought of exactly what you described @enki. Scale only at the beginning and first with simple multiplication and then do the rotation and translation chain. It also makes sense to transform the final versor to a 4x4 matrix before applying it to many points, since the computational effort is less.
I was still interested in the mathematical description of the outermorphism of scaling. Since I’m not a mathematician that took me a little while to get and what I’ll describe now may be trivial to you, but I think it might be helpful for someone like me who wants to understand it practically.
First, I think of scaling as moving whatever primitive you have (point, line, plane) from/to the origin by a factor. Since in PGA our 1-vectors are planes and we start with our outermorphism for the 1-vectors as @cgunn3 wrote, let’s have a look at a simple plane: x - 1 = 0, or E1 - E0 in PGA. This is the yz-plane moved by 1 along the x-axis. If we scale it by a factor s, it moves to x - s = 0, or E1 - s E0. Since this example is a moved yz-plane it doesn’t make a difference if you do uniform scaling or non-uniform scaling along x. Since any scaling of the 1-vector doesn’t change the plane, you could also write E1 / s - E0 which hints at how you can do non-uniform scaling.
We find that we can scale any plane a E0 + b E1 + c E2 + d E3 as follows: a E0 f + b E1 f / x + c E2 f / y + d E3 f / z. Here x, y and z are the non-uniform scaling factors and f is arbitrary as long as it is not zero. We will look at possible choices for f further below, but what I already want to mention here is that in general the resulting (multi-)vector will not be normalized after scaling even if it was before.
If you normalize after scaling though, f will be gone.
Now for the outermorphism, let’s see what happens to points, since that’s what I’m mostly interested in, so let’s scale the point a E021 + b E013 + c E032 + E123:
a E021 f³ / (x y) + b E013 f³ / (x z) + c E032 f³ / (y z) + E123 f³ / (x y z)
After normalization, we end up with the expected a E021 z + b E013 y + c E032 x + E123.
Now for the choice of f: you can of course simply choose 1, to save some computation if you implement this.
An alternative would be to set f = (x y z)^(1/3), since with this value you would get the scaled point immediately in the above example.
However, probably the smartest choice for f is f = x y z since this turns scaling the plane to a E0 x y z + b E1 y z + c E2 x z + d E3 x y. This avoids any divisions and thus also accounts for the case where one of the scaling factors is zero.
Lastly, I wanted to have a look at normal vectors, since these frequently cause problems when we try to implement shaders with non-uniform scaling. Representing normal vectors as ideal points (a E021 + b E013 + c E032) and scaling these doesn’t work. You have to take the dual view and consider them as planes (a E3 + b E2 + c E1) which can then be scaled and normalized.
I hope this helps someone.
Cheers