Hi Simon,
I have had experiences working when scripts to calculate least square solution start becoming computationally heavy and slow. You have done a good job on pinpointing the line taking the most amount of time.
The performance bottleneck in the calculation “err=psi'*(I-H*t)*psi” is due to the explicit formation of large intermediate matrices. If the observation matrix “H” has dimensions “m*n” and the measurement vector “psi” has length “m”, then “I” and “H*”t are both large” m*m” matrices. When “m” is large, the creation and multiplication of these matrices are computationally expensive and consume significant memory.
The expression calculates the sum of squared residuals. A more efficient computation can be achieved by reformulating the problem. The error, “err”, is the squared norm of the residual vector “r”, where “r” is the difference between the actual measurements “psi” and the values predicted by the model, “H.PCO”. This approach avoids large matrix-matrix operations. Additionally, the parameter estimate “PCO” can be calculated more efficiently and with better numerical stability. Instead of forming the normal equations with “t=(H'*H)\H' ”, one can use the MATLAB backslash operator (\) to solve the least-squares problem directly.
The following function implements these optimizations:
function [PCO,err] = calcPC_LSQ_optimized(H,psi)
% Obtain the least-squares solution directly using mldivide.
PCO = H\psi;
% Calculate the error as the sum of squared residuals.
residual = psi - H*PCO;
err = residual' * residual;
end
This updated method is more efficient as it relies on matrix-vector products and vector operations, which are significantly faster.
More details on solving linear least-squares problems can be found in the MATLAB documentation for - https://www.mathworks.com/help/matlab/ref/mldivide.html
