In[]:=
deploy
Fri 24 Jun 2022 11:37:40
In[]:=
cov[A_]:=A.A;cov[A_,B_]:=A.B;(*columnofcoefficientsneededtopredicti'thcoordinateofdata,datapointsarearrangedasrowsof"data"*)phi[i_,data_]:=Module[{Y1,Y2,p,residuals},Y1=Drop[data,{i}];(*independentvars,i'thcolumndropped*)Y2=data[[All,{i}]];(*targetvariable,i'thcolumn*)(*optimalpredictioncoeffs,3.26ofPourahmadi*)p=Inverse[Y1.Y1].(Y1.Y2);(*Formresiduals,checkthattheyareuncorrelatedwithoriginalpredictorvariables*)residuals=Y2-Y1.p;Assert[Total[cov[residuals,Y1],2]≐0];Insert[p,{0},{i}](*0coefficientinfrontofcurrentvariable*)];getB[X_]:=Table[phi[i,X]//Flatten,{i,n}];(*convertprecisionmatrixtomatrixofcoefficients*)invToPhi[inv_]:=Module[{D2,iiMat},D2=DiagonalMatrix@Diagonal@inv;IdentityMatrix[Length[inv]]-inv.Inverse[D2]];error[mat1_,mat2_]:=Print["full-rank"];X={{2,-2,2},{-1,1,1},{-2,-2,1},{1,1,3}};{m,n}=Dimensions[X];Print[getB[X]//MatrixForm]Print[invToPhi[Inverse[X.X]]//MatrixForm]Print["error: ",error[invToPhi@Inverse[X.X],getB[X]]];Print["rank-deficient + pseudo-inverse"];X={{2,-2,2},{-1,1,1},{-2,-2,1},{1,1,3}};{m,n}=Dimensions[X];Print[getB[X]//MatrixForm]Print[invToPhi[PseudoInverse[X.X]]//MatrixForm]Print["error: ",error[invToPhi[PseudoInverse[X.X]],getB[X]]//N]Print["rank-deficient + random noise"];X={{2,-2,2},{-1,1,1},{-2,-2,1},{1,1,3}};{m,n}=Dimensions[X];X2=X+0.00001*RandomVariate[NormalDistribution[],{m,n}];Print[getB[X]//MatrixForm]Print[invToPhi[Inverse[X2.X2]]//MatrixForm]Print["error: ",error[invToPhi[Inverse[X2.X2]],getB[X]]//N](*orthoprojected*)Print["orthorpojected"];xToPhi[X_]:=Module[{ii,P},ii=IdentityMatrix[Last@Dimensions@X];P=ii-PseudoInverse[X].X;(*orthoprojectorfornullspaceofX*)ii-P.Inverse[DiagonalMatrix[Diagonal[P]]]];Print["error orthoproject: ",error[xToPhi[X],getB[X]]//N]
2
Norm[mat1-mat2,"Frobenius"]
full-rank
0 | 4 67 | 2 5 |
4 73 | 0 | - 1 5 |
20 73 | - 10 67 | 0 |
0 | 4 67 | 2 5 |
4 73 | 0 | - 1 5 |
20 73 | - 10 67 | 0 |
error: 0
rank-deficient + pseudo-inverse
0 | - 1 2 | 7 4 | 7 8 |
-2 | 0 | 7 2 | 7 4 |
4 7 | 2 7 | 0 | - 1 2 |
8 7 | 4 7 | -2 | 0 |
0 | 1774 2237 | 53 2166 | 472 2833 |
887 1474 | 0 | - 1153 4332 | - 1786 2833 |
53 1474 | - 1153 2237 | 0 | 241 2833 |
118 737 | - 1786 2237 | 241 4332 | 0 |
error: 40.1094
rank-deficient + random noise
0 | - 1 2 | 7 4 | 7 8 |
-2 | 0 | 7 2 | 7 4 |
4 7 | 2 7 | 0 | - 1 2 |
8 7 | 4 7 | -2 | 0 |
0. | -0.499996 | 1.74995 | 0.874998 |
-2.00002 | 0. | 3.49994 | 1.75001 |
0.571444 | 0.285719 | 0. | -0.500012 |
1.14286 | 0.571425 | -1.99995 | 1.11022× -16 10 |
error: 9.65965×
-9
10
orthorpojected
error orthoproject: 0.
Orthoprojector method for zero columns
Orthoprojector method for zero columns
In[]:=
X={{2,-2,2,0,0},{-1,1,1,0,0},{-2,-2,1,0,0},{1,1,3,0,0}};X//MatrixForm
Out[]//MatrixForm=
2 | -2 | 2 | 0 | 0 |
-1 | 1 | 1 | 0 | 0 |
-2 | -2 | 1 | 0 | 0 |
1 | 1 | 3 | 0 | 0 |
In[]:=
P=IdentityMatrix[5]-PseudoInverse[X].X;P//MatrixForm
Out[]//MatrixForm=
0 | 0 | 0 | 0 | 0 |
0 | 0 | 0 | 0 | 0 |
0 | 0 | 0 | 0 | 0 |
0 | 0 | 0 | 1 | 0 |
0 | 0 | 0 | 0 | 1 |