This project is for educational purposes and can be done more directly by existing functionality. This post refers more to students trying to figure out how to solve the following problem for the first time
I'm lazier than you.
data = {{-1, 3}, {1, 1}, {2, 5}}; lm = LinearModelFit[data, {x, x^2}, x, IncludeConstantBasis -> False, WorkingPrecision -> \[Infinity]]; lm["BestFitParameters"] (* {-(12/11), 20/11} *)
Might want to mention built-in capabilities as well as the PseudoInverse function. Or at least note that this is for educational purposes and can be done more directly by existing functionality.
PseudoInverse