%% Total Least Squares (TLS): multivariate example with 5 random variables
% © 2022, Antonio Sala Piqueras, Universitat Politècnica de València. All rights reserved.
% This code successfully executed on Matlab R2022a
% Objective: illustrating a Total Least Squares case study involving five correlated variables.
% Presentations in video:
% http://personales.upv.es/asala/YT/V/tls51EN.html , http://personales.upv.es/asala/YT/V/tls52EN.html .
% Table of Contents
% Data generation (we are God)
% Ordinary LS (biased)
% TLS estimate
%     Data in rows
%     Data arranged in "columns"
% Addenda: what would happen if my scaling is wrong?... well, TLS won't work well!
% Conclusions
%% Data generation (we are God)
% Let us consider a "hidden" model where two variables (y_1,y_2) are a linear function of other three ones (x_1,x_2,x_3). Well, TLS will seek finding correlations in the vector (y_1,y_2,x_1,x_2,x_3).
% The model will be in the form
% (y_1\ y_2)_{N\times 2}=(x_1\ x_2\ x_3)_{N\times 3}\cdot \Theta_{3\times 2}
% so each sample of data has the variables in a "row", being N the number of samples.
% We may express an equivalent model in "column" form, of course, just transposing.

Xclean=randn(90000,3)*diag([25 1.7 18]);%random X data
ThetaClean=[2 40 5;-2 10 3]'
Yclean=Xclean*ThetaClean;
%TLS assumes clean linearly-related data are corrupted by independent
%measurement noise in all channels, let's do it...
stdx1=6; stdx2=0.4; stdx3=5; stdy1=0.5; stdy2=7.5;
stdX=diag([stdx1 stdx2 stdx3]); stdY=diag([stdy1 stdy2]);
X=Xclean+randn(size(Xclean))*stdX; %measurement noise in X
Y=Yclean+randn(size(Yclean))*stdY;%measurement noise in Y
%of course, only X and Y are actually "accessible" to the data engineer in the code below...

%
%% Ordinary LS (biased)
% Y=X\theta, we may obtain \theta with:

Th_LS_biased=pinv(X)*Y
ThetaClean

% Biased because we do not recover the correct theta even with "a lot" of data.
%
%% TLS estimate
% *We'll assume that standard deviations of the random additive noise are known... that assumption may actually be difficult in practice!.
%% Data in rows

size(X)
size(Y)
stdX=diag([stdx1 stdx2 stdx3]); %we repeat because we cannot access (supposedly) the data generation code.
stdY=diag([stdy1 stdy2]); %we repeat because we cannot access (supposedly) the data generation code.
Xesc=(X-mean(X))*inv(stdX); %we scale to zero mean and unit variance measurement noise (NOT unit variance of the data!)
Yesc=(Y-mean(Y))*inv(stdY); %we scale to zero mean and unit variance measurement noise (NOT unit variance of the data!)
Data_esc=[Yesc Xesc];  %Now all five columns have "unit variance measurement noise"
[N,m]=size(Data_esc)

% *Note: in the video, I first divided by standard deviation and I subtracted the mean later on... It's OK but, well, what everybody does to center data is first subtracting the mean and then dividing by standard deviation, so I changed it in the materials.

tic
[U,S,V]=svd(Data_esc/sqrt(N-1),'econ'); %TLS and SVD are the same with the proposed scaling.
%dividing by sqrt(N-1), S has units of standard deviation
toc %svd is fast.

% The standard deviation of each principal component is:

diag(S)' %I should find the model as "unit noise" singular values, given the above scaling
size(U)
size(V) %the 5x5 matrix is the one to use in TLS
Model_scaled=V(:,4:5) %model relating scaled variables

% So we have "[Yesc Xesc]*Model_scaled \approx 0". TLS does not know whether any of these variables are "inputs" or "outputs". As we generated data with Y=X\theta, let us solve for Y to check what we get. Once we undo the scaling, we will have
% YM_1 + XM_2=0,    thus      YM_1=-XM_2, and finally Y=-X\cdot M_2M_1^{-1}.

ModX_sc=Model_scaled(3:5,:) 
ModY_sc=Model_scaled(1:2,:)  
ModX=stdX\ModX_sc %matrix multiplying X in original units in the found model (M2 above)
ModY=stdY\ModY_sc  %matrix multiplying Y in original units in the found model (M1 above)
ThetaEstimTLS=-ModX*inv(ModY)  %estimated TLS parameter, solving for Y explicitly

% It's unbiased when comparing to our ideal "clean" parameter we used when generating the data:

ThetaClean

%
%%
%% Data arranged in "columns"
%

X=X';Y=Y'; 
size(X)
size(Y)

% We would now have

stdX=diag([stdx1 stdx2 stdx3]);
stdY=diag([stdy1 stdy2]);
Xesc=inv(stdX)*X;
Yesc=inv(stdY)*Y;
Data_esc=[Yesc; Xesc]; 
[m,N]=size(Data_esc)
Data_esc=Data_esc-sum(Data_esc,2)/N; %subtract the mean, now summing along dimension 2
[U,S,V]=svd(Data_esc/sqrt(N-1),'econ'); 
diag(S)'
size(U) %choose 5x5 for TLS
Model_scaled=U(:,4:5)' %new model, now it's the transpose of the earlier one, but it's the same.
ModX_sc=Model_scaled(:,3:5);
ModY_sc=Model_scaled(:,1:2);  
ModX=ModX_sc/stdX  
ModY=ModY_sc/stdY 
ThetaEstimTLS=-inv(ModY)*ModX %this is how we now must solve Y=\theta X
ThetaClean'

% So it's irrelevant whether I arrange data in rows or columns, results are equivalent, of course.
%
%% Addenda: what would happen if my scaling is wrong?... well, TLS won't work well!

DataWrong=[Y; X]; %NO scaling, because I (wrongly) thought measurement noise was already "equalised".
DataWrong=DataWrong-sum(DataWrong,2)/N;
[U,S,V]=svd(DataWrong/sqrt(N-1),'econ');
diag(S)' %only one equation relating things? (significantly smaller singular value)?
Model_wrong_scaling=U(:,4:5)' %no need of undoing the scaling step, because there was no such step: original units already
Theta_wrong_scaling=-inv(Model_wrong_scaling(:,1:2))*Model_wrong_scaling(:,3:5) %biased
ThetaClean'
idealModel=[eye(2) -ThetaClean'];
Model_good_scaling=[ModY ModX] %the "good TLS one", in original un-scaled coordinates

% Saying that Model*Data\approx 0 amounts to saying that the allowed data according to the model are the subspace given by null(Model). Let us compute the "angle" between these subspaces to summarise the modelling error in a single number:

subspace(null(Model_wrong_scaling),null(idealModel))*180/pi %angle in degrees
subspace(null(Model_good_scaling),null(idealModel))*180/pi %angle in degrees
subspace(null([-eye(2) Th_LS_biased']),null(idealModel))*180/pi

% So, the "angle" bias due to wrong TLS scaling is even worse than the bias of the standard LS solution. Wow! Nevertheless, the bias less than 3 degrees might not be so much in applications (I would have, say, more severe "variance-related" problems if I had 50 data instead of 90000). Note however, that "subspace angle" depends on the "chosen" scaling, so changing a variable from volts to milivolts would change the angle output above.
%% Conclusions
% If data are "very abundant", and generated in accordance with theory and I know the "measuring noise variance", then everything works as theory predicts.
% If I have not so many data, or I make a wrong guess with the scaling, then TLS may be biased. Real life is hard.
%
%