Yes. I believe that as far as it goes, that specific statement is correct. But you are right that I was neglecting the truncation in my example. I think it is harder to show backward instability impacting forward errors with a simple example than I thought. That’s a little embarrassing.
The fact that explicit (pseudo)inverses are not backward stable is something that applies even in the special case of a square invertible matrix. There’s a nice discussion of this on page 260 of Higham’s Accuracy and Stability of Numerical Algorithms (2nd edition). I don’t have a nice pdf or scan to include, but a typed in quote from that is "Suppose X = A^{-1} is formed exactly, and that the only rounding errors are in forming x = \mbox{fl}(Xb). Then \hat{x} = (X+\Delta X)b where |\Delta X| \leq \gamma_n |X|, by (3.11). So A\hat{x} = A(X+\Delta X)b = (I+A\Delta X) b and the best possible residual bound is
Since…" Here \gamma_n = nu /(1-nu) where u is the unit roundoff.
It goes on to note that that residual bound is not consistent with backward stability if \|x\|_\infty \ll \||A^{-1}| |b|\|_\infty. So even getting an exact inverse and multiplying by b with numerical errors is not going to be backward stable. This is the source of potentially large residuals in ill-conditioned square systems solved using an inverse.
For the section from Trefethen and Bau that you quoted, moving the parenthesis results in an explicit computation of an inverse and this sort of argument applies in the square invertible case. Even forming V*(S\U') exactly isn’t good enough for backward stability once rounding in the matrix vector multiply (represented by \Delta X) is taken into account.
A square system is obviously a simpler special case of the pseudoinverse. For a full rank (but ill conditioned) least squares problem I would have thought that you would be able to observe condition squaring without too much effort. That’s what I was trying for in that example. My reasoning was that the conditioning for a LS problem involves a term proportional to both \kappa^2(A) and \|r\|. So the condition squaring does not impact small residual problems so much if you use a backward stable algorithm. This is the forward accuracy advantage of factorization over the normal equations. The pseudoinverse would not be computed using the normal equations, but the columns of the pseudoinverse are solutions to the LS problem with b=e_k. Even if the given problem has a small residual that damps out condition squaring, the columns of the pseudoinverse will not necessarily correspond to problems with a small residual. So I would expect errors proportional to \kappa^2(A)u to be visible in a solution formed from a linear combination of the columns of the pseudoinverse without any mitigation due to a small residual. In a backward stable solution those error terms would be multiplied by \|r\|. Since I forgot about the automatic truncation, that’s what I thought I was observing. I didn’t look at the results as closely as I should have once they seemed to confirm my expectations.
Possibly there’s a flaw in that (admittedly loose) reasoning or it’s simply harder to come up with a nice example illustrating the effect of instability on forward errors than I thought. I’m not sure which now…