It depends a lot of the problem. For a general least squares solution of the Ax=b problem, is your A matrix dense, sparse, structured? Do you need to solve it very accurately or just approximately? Are you going to be using this over and over with the same A but different b's? In that case, you may profit from a one time Cholesky factorization that can be used over and over very quickly. Is your linear operator A explicit (as in you can actually inspect the matrix) or is it a linear operator function in which you can only see the result of A(x) (and its adjoint A_adj(y)) ? In this case, you cannot use common factorizations (QR, cholesky, svd) but will have to use iterative methods (ie. Krylov space methods).

There's a whole branch in applied math that's called numerical linear algebra and deals with all these kinds of situations. One of the most important results is that you should almost never do an explicit computation of the inverse of a matrix. It is expensive and can be highly unstable. If you ever see a A\b or pinv(A)*b in somebody's code, it should raise a big red flag.

Instead of

  𝛉 = pinv(𝐗)*y

write

  𝛉 = y\𝐗

It will automatically dispatch to the most appropriate algorithm depending on the type of X (Cholesky, QR, LU/Gauss, …)

Often one is not really interested in solving Ax=b for all possible b but just one particular b. In that case one can iteratively solve for x by algorithms like conjugate gradient which only require multiplying A times the current best guess for x, which can be very fast and avoids numerical instability issues in many cases.

QR, Cholesky

QR and Cholesky factorization (QR is the most versatile and is always stable, so usually a better first choice)