Project : Propogation of Errors in Decompositions and Distributions of Errors/Q’s

Project Reports

Project Report #1 (3/23):

So far our group has looked into the Error distribution of the QR as well as the distributions of the elements of Q, the elements of Q were found to be non - normal and we' re looking into literature to see if anyone has found out what exactly they are.We' ve also developed code for the QR and GS decompositions (LU is in the works) which degrade the accuracy during the matrix decomposition and are writing code to do a large simulation to view the effects of these errors upon the decomposition.

Project Report #2 (3/28):

Investigation of the q’s in the QR Decomposition and their Distribution:

To investiage the distribution of the entries of the Q in the QR decomposition a matrix of size 100 x 100 was generated and a normal distribution was fitted to the entries with a histogram showing them overlaying and a a quantile plot of the values as well as a Chi-Square Test for normality. As you can see it looks quite plausibly normal...

The P-value of the Chi-Square test is .918688, which gives no cause to doubt that the distribution is non-normal, however it the hypothesis that the Distribution is normal is true, then the distribution of the P-values should be uniform. So the previous test was repeated 3000 times and we check whether the distribution of the P’s looks uniform...

n=100;
reps=3000;
Timing[pvalues=Table[A=GenMatrixConditionedAs[{n,n},10];{Q,R}=QR[A];PearsonChiSquareTest[Flatten[Q]],{i,reps}];];
Histogram[pvalues,Automatic,"PDF",ImageSize->size,PlotLabel->"Histogram of P-values of Chi-square tests"]
Print["P-value of Chi-Square Test checking p-values look Uniform: p = ",PearsonChiSquareTest[pvalues,UniformDistribution[{0,1}]]];

So the Distribution of those P-values is highly unlikely to be uniform and we can conclude that in general the enteries of the Q are not normally distributed. However as you can recall, then columns treated as samples look more and more normal as you go from left to right, this is going to be explored later, though it is still very unlikely that the final column is normally distributed as the values cannot be independent of one another due to the fact that the norm of the column must be 1.

Investigation of the Patterns in the errors of the QR and GS decompositions

We’ve just realized that our GS code isn’t double orthogonalizing, but we intend to fix it. In the mean time you’ll find below our investigations of what happens when there is a single entry in R as you construct the QR decomposition as well as when you miscast the Orthogonalization of the GS code down to a single from a double. 

It appears that for the QR decomposition, when you want to solve R, then answer is only changed after the point at which the error was made in the matrix.

And for the GS, the errors accumulate as the algorithm works its way through A, which is expected, and that the enteries in the resultant Q are not normally distributed.

We’re still looking into what exactly the distribution of the Errrors of the decomposition actually is, the raw errors look like the Laplace Distribution (double exponential), but the Pearson Chi-square test gives an exceedingly small p-value which leads one to conclude that though it looks that way, it actually isn’t.

You can see the work we’ve done below.
↓

Functions used in the document

Code for degrading the accuracy of a double

(* This Implements the Truncation Rounded Scheme *)
Degrade[number_]:=Module[{n},
n=RealDigits[number,2];
n={n[[1, 1 ;; 23]],n[[2]]};
Sign[number] N[FromDigits[n,2]]
];

(* This Implements the "Round to Nearest, Tie to Even" Rounding Scheme *)
Degrade2[number_]:=Module[{n,chop=23,i,m},
If[number≠0,
m=RealDigits[number,2]; 
n=m[[1]];
If[n[[chop+1]]==1 && Total[n[[chop+2;;-1]]]≥1,
	n[[chop]]=n[[chop]]+1;
	Do[
		If[n[[i+1]]==2,
			n[[i]]=n[[i]]+1;
			n[[i+1]]=0;
		];
	,{i,chop-1,1,-1}];
	If[n[[1]]==2,n[[1]]=1];
];

n=Sign[number] N[FromDigits[{n[[1;;chop]],m[[2]]},2]],
n=number;
];
n
];

(* This Implements the "Round to Nearest, Tie to Odd" Rounding Scheme *)
Degrade3[number_]:=Module[{n,chop=23,i,m},
m=RealDigits[number,2]; 
n=m[[1]];
If[n[[chop+1]]==1,
	n[[chop]]=n[[chop]]+1;
	Do[
		If[n[[i+1]]==2,
			n[[i]]=n[[i]]+1;
			n[[i+1]]=0;
		];
	,{i,chop-1,1,-1}];
	If[n[[1]]==2,n[[1]]=1];
];

Sign[number] N[FromDigits[{n[[1;;chop]],m[[2]]},2]]
];

(* Degrades all of the enteries in a Vector or Matrix (Any List object) using Degrade2 *)
DegradeM[A_List]:=Module[{},
Map[Degrade2,A,{-1}]
];

x=-N[.2];
FullForm[{Degrade[x],x}]
Plot[((Degrade[x]-x)/x),{x,0.1,1}]
FullForm[{Degrade2[x],x}]
Plot[((Degrade2[x]-x)/x),{x,0.1,1}]
FullForm[{Degrade3[x],x}]
Plot[((Degrade3[x]-x)/x),{x,0.1,1}]

Code for Generating Badly Conditioned Matricies

(* FIND THE CONDITION NUMBER OF A MATRIX *)
ConditionNumber[A_]:=Module[{w},
w=SingularValueList[A];
Max[w]/Min[w]
]

(* GENERATING AN UNSPECIFIED MATRIX OF SPECIFIED CONDITION AND STRUCTURE *)
GenMatrixConditionedAs[{n0_,m0_},conditionNumber0_:1,symmetricCheck0_:0,positiveDefCheck0_:0]:=Module[{n=n0,m=m0,conditionNumber=conditionNumber0,symmetricCheck=symmetricCheck0,positiveDefCheck=positiveDefCheck0,A,u,v,w,diagLength=Min[{n0,m0}],e},
(* Random Matrix Appropriate Size *)
A=RandomReal[{-1,1},{n,m}];
(* Take the SVD of the random Matrix to have bases for the col/row spaces *)
{u,w,v}=SingularValueDecomposition[A];
(* Constructing the Singular/Eigen Values of the new Matrix *)
e=Table[
negative=RandomInteger[]; (* Choosing whether it's positive or negative *)
(* Make sure it's positive if flagged for symmetric pos. def. *)
If[positiveDefCheck==1&&symmetricCheck==1&&n==m,negative=0]; 
(* Equally spaced singular/eigen values from 1 to condition number *)
(1+((conditionNumber-1)/(diagLength-1))i)*(((-1))^negative),{i,0,diagLength-1}
];
(* Put the values in the center matrix W *)
w=SparseArray[Band[{1,1}]->e,{n,m}];
(* Checking if flagged for symmetry *)
If[n==m && symmetricCheck==1,
(* If Symmetric *)
A=u.w.Transpose[u];,
(* If not symmetric *)
A=u.w.ConjugateTranspose[v];
];
A
];


GenMatrixConditionedAsN[{n0_,m0_},conditionNumber0_:1,symmetricCheck0_:0,positiveDefCheck0_:0]:=Module[{n=n0,m=m0,conditionNumber=conditionNumber0,symmetricCheck=symmetricCheck0,positiveDefCheck=positiveDefCheck0,A,u,v,w,diagLength=Min[{n0,m0}],e},
(* Random Matrix Appropriate Size *)
A=RandomReal[NormalDistribution[],{n,m}];
(* Take the SVD of the random Matrix to have bases for the col/row spaces *)
{u,w,v}=SingularValueDecomposition[A];
(* Constructing the Singular/Eigen Values of the new Matrix *)
e=Table[
negative=RandomInteger[]; (* Choosing whether it's positive or negative *)
(* Make sure it's positive if flagged for symmetric pos. def. *)
If[positiveDefCheck==1&&symmetricCheck==1&&n==m,negative=0]; 
(* Equally spaced singular/eigen values from 1 to condition number *)
(1+((conditionNumber-1)/(diagLength-1))i)*(((-1))^negative),{i,0,diagLength-1}
];
(* Put the values in the center matrix W *)
w=SparseArray[Band[{1,1}]->e,{n,m}];
(* Checking if flagged for symmetry *)
If[n==m && symmetricCheck==1,
(* If Symmetric *)
A=u.w.Transpose[u];,
(* If not symmetric *)
A=u.w.ConjugateTranspose[v];
];
A
];

Code for a QR Decomposition

Householder[x_List,i_]:= Module[{n=Length[x],u,H, τ=Norm[x[[i;;-1]]]},
H=SparseArray[Band[{1,1}]->1,{n,n}];
τ=Sign[x[[i]]] τ; 
u=Join[Table[0,{i-1}],{1},((x[[i+1;;-1]])/(τ+x[[i]]))];
H-(2/((Norm[u])^2))Outer[Times,u,u]
];
HouseholderDense[x_List,i_]:= Module[{n=Length[x],u,H, τ=Norm[x[[i;;-1]]]},
H=IdentityMatrix[n];
τ=Sign[x[[i]]] τ; 
u=Join[Table[0,{i-1}],{1},((x[[i+1;;-1]])/(τ+x[[i]]))];
H-(2/((Norm[u])^2))Outer[Times,u,u]
];

QR[A_List]:=Module[{R=A,n,m,Q,Qtemp,k},
{n,m}=Dimensions[A];
Q=IdentityMatrix[n];
If[n>m,k=m,k=n-1];
Do[
Qtemp=Householder[R[[All,j]],j];
R=Qtemp.R;
Q=Q.Qtemp;
,{j,1,k}];
{Q,R}
];

Code for a QR Decomposition with loss of precision in R at (50,50) when k=50 in the middle of the QR

QRloss1[A_List]:=Module[{R=A,n,m,Q,Qtemp,k},
{n,m}=Dimensions[A];
Q=IdentityMatrix[n];
If[n>m,k=m,k=n-1];
Do[
If[j==50, R[[50,50]]=Degrade2[R[[50,50]]]];
Qtemp=Householder[R[[All,j]],j];
R=Qtemp.R;
Q=Q.Qtemp;
,{j,1,k}];
{Q,R}
];

QRloss2[A_List]:=Module[{R=A,n,m,Q,Qtemp,k},
{n,m}=Dimensions[A];
Q=IdentityMatrix[n];
If[n>m,k=m,k=n-1];
Do[
If[j==50, R[[60,50]]=Degrade2[R[[60,50]]]];
Qtemp= Householder[R[[All,j]],j];
R=Qtemp.R;
Q=Q.Qtemp;
,{j,1,k}];
{Q,R}
];

QRloss3[A_List]:=Module[{R=A,n,m,Q,Qtemp,k},
{n,m}=Dimensions[A];
Q=IdentityMatrix[n];
If[n>m,k=m,k=n-1];
Do[
Qtemp=Householder[R[[All,j]],j];
If[j==50, Qtemp[[50,50]]=Degrade2[Qtemp[[50,50]]]];
R=Qtemp.R;
Q=Q.Qtemp;
,{j,1,k}];
{Q,R}
];

QRloss4[A_List]:=Module[{R=A,n,m,Q,Qtemp,k},
{n,m}=Dimensions[A];
Q=IdentityMatrix[n];
If[n>m,k=m,k=n-1];
Do[
Qtemp=Householder[R[[All,j]],j];
R=Qtemp.R;
If[j==50, Qtemp[[50,50]]=Degrade2[Qtemp[[50,50]]]];
Q=Q.Qtemp;
,{j,1,k}];
{Q,R}
];

Code for a LU Decomposition

LU[A_]:=Module[{LU=A, n=Length[A],d,e,r,i,j},
Do[
d=1/LU[[r,r]];
Do[
e = d LU[[i,r]];
LU[[i,r]]=e;
Do[
LU[[i,j]]-= e* LU[[r,j]],
{j,r+1,n}
],
{i,r+1,n,1}],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
];

LUPivot[A_]:=Module[{LU=A, n=Length[A],perm,pivot,d,e,r,i,j},
perm=Table[i,{i,n}];
Do[
(* Partial Pivoting *)
pivot=Ordering[Abs[LU[[r;;n,r]]],-1][[1]]+r-1;
perm[[{r,pivot}]]=perm[[{pivot,r}]];
LU[[{r,pivot}]]=LU[[{pivot,r}]];
(* End Partial Pivoting *)
d=1/LU[[r,r]];
Do[
e = d LU[[i,r]];
LU[[i,r]]=e;
Do[
LU[[i,j]]-= e* LU[[r,j]],
{j,r+1,n}
],
{i,r+1,n,1}],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU], perm}
]

Code for GS Decomposition

(*Gram-Schmidt with Double orthogonalization*)
GSO[A_]:=
Module[
{T=A,r,rHat,k,i,j,m=Length[(A^\[Transpose])]},
r=RandomReal[{-1,1},{m,m}]*0;(*Initializing r*)
For[k=1,k≤m,k++,
For[i=1,i≤k-1,i++,
r[[i,k]]=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*r[[i,k]](*first orthogonalization*);
];
For[i=1,i≤k-1,i++,
rHat=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*rHat(*second orthogonalization*);
r[[i,k]]+=rHat;
];
r[[k,k]]=Norm[T[[All,k]]];
If[r[[k,k]]≤$MachineEpsilon,
Return["The Matrix does not have a full rank!"]];
T[[All,k]]=(1/r[[k,k]])*T[[All,k]];
];
{T,r}(*T=Q, r=R*)
];

Code for GS Decomposition without normalization on the last step

(*Gram-Schmidt with Double orthogonalization*)
GSOn[A_]:=
Module[
{T=A,r,rHat,k,i,j,m=Length[(A^\[Transpose])]},
r=RandomReal[{-1,1},{m,m}]*0;(*Initializing r*)
For[k=1,k≤m,k++,
For[i=1,i≤k-1,i++,
r[[i,k]]=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*r[[i,k]](*first orthogonalization*);
];
For[i=1,i≤k-1,i++,
rHat=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*rHat(*second orthogonalization*);
r[[i,k]]+=rHat;
];
r[[k,k]]=Norm[T[[All,k]]];
If[r[[k,k]]≤$MachineEpsilon,
Return["The Matrix does not have a full rank!"]];
If[k≠m,
T[[All,k]]=(1/r[[k,k]])*T[[All,k]];
];
];
{T,r}(*T=Q, r=R*)
]

Code for QR with missed cast

Householder[x_List,i_]:= Module[{n=Length[x],u,H, τ=Norm[x[[i;;-1]]]},
H=SparseArray[Band[{1,1}]->1,{n,n}];
τ=Sign[x[[i]]] τ; 
u=Join[Table[0,{i-1}],{1},((x[[i+1;;-1]])/(τ+x[[i]]))];
H-(2/((Norm[u])^2))Outer[Times,u,u]
];

Householder2[x_List,i_]:= Module[{n=Length[x],u,H, τ=Norm[x[[i;;-1]]]},
H=SparseArray[Band[{1,1}]->1,{n,n}];
τ=Sign[x[[i]]] τ; 
u=Join[Table[0,{i-1}],{1},((x[[i+1;;-1]])/(τ+x[[i]]))];
H-DegradeM[(2/((Norm[u])^2))Outer[Times,u,u]]
];


QRloss1[A_List]:=Module[{R=A,n,m,Q,Qtemp,k},
{n,m}=Dimensions[A];
Q=IdentityMatrix[n];
If[n>m,k=m,k=n-1];
Do[
Qtemp=DegradeM[Householder[R[[All,j]],j]];
R=Qtemp.R;
Q=Q.Qtemp;
,{j,1,k}];
{Q,R}
];
QRloss2[A_List]:=Module[{R=A,n,m,Q,Qtemp,k},
{n,m}=Dimensions[A];
Q=IdentityMatrix[n];
If[n>m,k=m,k=n-1];
Do[
Qtemp=Householder[R[[All,j]],j];
R=DegradeM[Qtemp.R];
Q=DegradeM[Q.Qtemp];
,{j,1,k}];
{Q,R}
];

QRloss3[A_List]:=Module[{R=A,n,m,Q,Qtemp,k},
{n,m}=Dimensions[A];
Q=IdentityMatrix[n];
If[n>m,k=m,k=n-1];
Do[
Qtemp=Householder2[R[[All,j]],j];
R=Qtemp.R;
Q=Q.Qtemp;
,{j,1,k}];
{Q,R}
];

Code for LU with missed cast

LUloss1[A_]:=Module[{LU=A, n=Length[A],d,e,r,i,j},
Do[
If[r==50,LU[[r,r]]=Degrade2[LU[[50,50]]];];
d=1/LU[[r,r]];
Do[
e = d LU[[i,r]];
LU[[i,r]]=e;
Do[
LU[[i,j]]-= e* LU[[r,j]],
{j,r+1,n}
],
{i,r+1,n,1}],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
]

LUloss2[A_]:=Module[{LU=A, n=Length[A],d,e,r,i,j,loss},
loss=0;
Do[
d=1/LU[[r,r]];
If[r==50&&loss==0,LU[[50,60]]=Degrade2[LU[[50,60]]];loss=1;];
Do[
e = d LU[[i,r]];
LU[[i,r]]=e;
Do[
LU[[i,j]]-= e* LU[[r,j]],
{j,r+1,n}
],
{i,r+1,n,1}],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
]

LUloss3[A_]:=Module[{LU=A, n=Length[A],d,e,r,i,j,loss},
loss=0;
Do[
d=1/LU[[r,r]];
Do[
e = d LU[[i,r]];
LU[[i,r]]=e;
If[r==50&&i==55&&loss==0,LU[[55,50]]=Degrade2[LU[[55,50]]]; loss=1;];
Do[
LU[[i,j]]-= e* LU[[r,j]],
{j,r+1,n}
];,
{i,r+1,n,1}],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
]

LUloss4[A_]:=Module[{LU=A, n=Length[A],d,e,r,i,j,loss},
loss=0;
Do[
d=1/LU[[r,r]];
Do[
e = d LU[[i,r]];
LU[[i,r]]=e;
Do[
LU[[i,j]]-= e* LU[[r,j]],
{j,r+1,n}
],
If[r==50&&i==55&&loss==0,LU[[55,50]]=Degrade2[LU[[55,50]]]; loss=1;];
{i,r+1,n,1}],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
]


LUloss5[A_]:=Module[{LU=A, n=Length[A],d,e,r,i,j,loss},
loss=0;
Do[
d=1/LU[[r,r]];
Do[
e = d LU[[i,r]];
LU[[i,r]]=e;
If[r==50&&i==50&&loss==0,LU[[50,45]]=Degrade2[LU[[55,50]]]; loss=1;];
Do[
LU[[i,j]]-= e* LU[[r,j]];
,{j,r+1,n}
],
{i,r+1,n,1}],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
]

LUPivotloss[A_]:=Module[{LU=A, n=Length[A],perm,pivot,d,e,r,i,j},
perm=Table[i,{i,n}];
Do[
(* Partial Pivoting *)
pivot=Ordering[Abs[LU[[r;;n,r]]],-1][[1]]+r-1;
perm[[{r,pivot}]]=perm[[{pivot,r}]];
LU[[{r,pivot}]]=LU[[{pivot,r}]];
(* End Partial Pivoting *)
d=1/LU[[r,r]];
Do[
e = d LU[[i,r]];
LU[[i,r]]=e;
Do[
LU[[i,j]]-= e* LU[[r,j]],
{j,r+1,n}
],
{i,r+1,n,1}],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU], perm}
]

Code for Cholesky Decomposition

OutChol[A_]:=Module[{dim=Dimensions[A],R=0 A,n,m,r=A[[1,1]],b=A[[1,2;;-1]],Ah=A[[2;;-1,2;;-1]]},
{n,m}=dim;
If[n==1,
R=(Sqrt[r]);
,
r=(Sqrt[r]);
b=b/r;
Ah=Ah-Outer[Times,b,b];
R[[1,1]]=r;
R[[1,2;;-1]]=b;
R[[2;;-1,2;;-1]]=OutChol[Ah];
];
R
];

Chol[A_]:=Module[{n=Dimensions[A][[1]]},
Block[{$RecursionLimit=n*4+1},
OutChol[A]]
];

Code for Cholesky Decomposition with missed cast

OutCholloss2[A_]:=Module[{dim=Dimensions[A],R=0 A,n,m,r=A[[1,1]],b=A[[1,2;;-1]],Ah=A[[2;;-1,2;;-1]]},
{n,m}=dim;
If[n==1,
R=(Sqrt[r]);
,
If[n==50,
b[[20]]=Degrade2[b[[20]]];
];
r=(Sqrt[r]);
b=b/r;
Ah=Ah-Outer[Times,b,b];
R[[1,1]]=r;
R[[1,2;;-1]]=b;
R[[2;;-1,2;;-1]]=OutCholloss1[Ah];
];
R
];

OutCholloss1[A_]:=Module[{dim=Dimensions[A],R=0 A,n,m,r=A[[1,1]],b=A[[1,2;;-1]],Ah=A[[2;;-1,2;;-1]]},
{n,m}=dim;
If[n==1,
R=(Sqrt[r]);
,
If[n==50,
r=Degrade2[r];
];
r=(Sqrt[r]);
b=b/r;
Ah=Ah-Outer[Times,b,b];
R[[1,1]]=r;
R[[1,2;;-1]]=b;
R[[2;;-1,2;;-1]]=OutCholloss2[Ah];
];
R
];

Cholloss[Ain_,k_]:=Module[{n=Length[Ain],A=Ain},
Block[{$RecursionLimit=n*4+1},
If[k==1,
A=OutCholloss1[A]];
If[k==2,
A=OutCholloss2[A]];
A
]
];

Code for a GS with a missed cast.

(*Code for QR decomposition with loss of precision in R at (25,25) when k=25 at the middle of the QR*)
GSO1[A_]:=
Module[
{T=A,r,rHat,k,i,j,m=Length[(A^\[Transpose])]},
r=RandomReal[{-1,1},{m,m}]*0;(*Initializing r*)
For[k=1,k≤m,k++,
For[i=1,i≤k-1,i++,
(*If[i==25,T[[25,25]]=Degrade2[T[[25,25]]]];*)
If[k==50,T[[50,50]]=Degrade2[T[[50,50]]]];
r[[i,k]]=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*r[[i,k]](*first orthogonalization*)
(*If[{i,k}=={25,25},T[[i,k]]=Degrade2[T[[i,k]]]]*)
];
For[i=1,i≤k-1,i++,
rHat=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*rHat(*second orthogonalization*);
r[[i,k]]+=rHat;
];
r[[k,k]]=Norm[T[[All,k]]];
If[r[[k,k]]≤$MachineEpsilon,
Return["The Matrix does not have a full rank!"]];
T[[All,k]]=(1/r[[k,k]])*T[[All,k]];
];
{T,r}(*T=Q, r=R*)
];

(*Code for GSO with missed cast*)
GSO2[A_]:=
Module[
{T=A,r,rHat,k,i,j,m=Length[(A^\[Transpose])],temp},
r=RandomReal[{-1,1},{m,m}]*0;(*Initializing r*)
For[k=1,k≤m,k++,
For[i=1,i≤k-1,i++,
temp=T[[All,k]].T[[All,i]];
If[k==50 && i==25,temp=Degrade2[temp];];
r[[i,k]]=temp;
T[[All,k]]-=T[[All,i]]*r[[i,k]](*first orthogonalization*);
];
For[i=1,i≤k-1,i++,
rHat=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*rHat(*second orthogonalization*);
r[[i,k]]+=rHat;
];
r[[k,k]]=Norm[T[[All,k]]];
If[r[[k,k]]≤$MachineEpsilon,
Return["The Matrix does not have a full rank!"]];
T[[All,k]]=(1/r[[k,k]])*T[[All,k]];
];
{T,r}(*T=Q, r=R*)
];

GSO3[A_]:=
Module[
{T=A,r,rHat,k,i,j,m=Length[(A^\[Transpose])]},
r=RandomReal[{-1,1},{m,m}]*0;(*Initializing r*)
For[k=1,k≤m,k++,
For[i=1,i≤k-1,i++,
r[[i,k]]=T[[All,k]].T[[All,i]];
If[k==50,r[[i,k]]=Degrade2[r[[i,k]]];];
T[[All,k]]-=T[[All,i]]*r[[i,k]](*first orthogonalization*);
];
For[i=1,i≤k-1,i++,
rHat=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*rHat(*second orthogonalization*);
r[[i,k]]+=rHat;
];
r[[k,k]]=Norm[T[[All,k]]];
If[r[[k,k]]≤$MachineEpsilon,
Return["The Matrix does not have a full rank!"]];
T[[All,k]]=(1/r[[k,k]])*T[[All,k]];
];
{T,r}(*T=Q, r=R*)
]
GSO4[A_]:=
Module[
{T=A,r,rHat,k,i,j,m=Length[(A^\[Transpose])]},
r=RandomReal[{-1,1},{m,m}]*0;(*Initializing r*)
For[k=1,k≤m,k++,
For[i=1,i≤k-1,i++,
r[[i,k]]=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*r[[i,k]](*first orthogonalization*);
];
For[i=1,i≤k-1,i++,
rHat=T[[All,k]].T[[All,i]];
If[k==50 && i==25,rHat=Degrade2[rHat];];
T[[All,k]]-=T[[All,i]]*rHat(*second orthogonalization*);
r[[i,k]]+=rHat;
];
r[[k,k]]=Norm[T[[All,k]]];
If[r[[k,k]]≤$MachineEpsilon,
Return["The Matrix does not have a full rank!"]];
T[[All,k]]=(1/r[[k,k]])*T[[All,k]];
];
{T,r}(*T=Q, r=R*)
]
GSO5[A_]:=
Module[
{T=A,r,rHat,k,i,j,m=Length[(A^\[Transpose])]},
r=RandomReal[{-1,1},{m,m}]*0;(*Initializing r*)
For[k=1,k≤m,k++,
For[i=1,i≤k-1,i++,
r[[i,k]]=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*r[[i,k]](*first orthogonalization*);
];
For[i=1,i≤k-1,i++,
rHat=T[[All,k]].T[[All,i]];
If[k==50, rHat=Degrade2[rHat];];
T[[All,k]]-=T[[All,i]]*rHat(*second orthogonalization*);
r[[i,k]]+=rHat;
];
r[[k,k]]=Norm[T[[All,k]]];
If[r[[k,k]]≤$MachineEpsilon,
Return["The Matrix does not have a full rank!"]];
T[[All,k]]=(1/r[[k,k]])*T[[All,k]];
];
{T,r}(*T=Q, r=R*)
]
GSO6[A_]:=
Module[
{T=A,r,rHat,k,i,j,m=Length[(A^\[Transpose])],tmp},
r=RandomReal[{-1,1},{m,m}]*0;(*Initializing r*)
For[k=1,k≤m,k++,
For[i=1,i≤k-1,i++,
r[[i,k]]=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*r[[i,k]](*first orthogonalization*);
];
For[i=1,i≤k-1,i++,
rHat=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*rHat(*second orthogonalization*);
r[[i,k]]+=rHat;
];
r[[k,k]]=Norm[T[[All,k]]];
If[r[[k,k]]≤$MachineEpsilon,
Return["The Matrix does not have a full rank!"]];
tmp=(1/r[[k,k]])*T[[All,k]];
If[k==50,tmp[[k]]=Degrade2[tmp[[k]]];];
T[[All,k]]=tmp;
];
{T,r}(*T=Q, r=R*)
];



GSO7[A_]:=
Module[
{T=A,r,rHat,k,i,j,m=Length[(A^\[Transpose])],tmp},
r=RandomReal[{-1,1},{m,m}]*0;(*Initializing r*)
For[k=1,k≤m,k++,
For[i=1,i≤k-1,i++,
r[[i,k]]=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*r[[i,k]](*first orthogonalization*);
];
For[i=1,i≤k-1,i++,
rHat=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*rHat(*second orthogonalization*);
r[[i,k]]+=rHat;
];
If[k==50,T[[k,k]]=Degrade2[T[[k,k]]];];
r[[k,k]]=Norm[T[[All,k]]];
If[r[[k,k]]≤$MachineEpsilon,
Return["The Matrix does not have a full rank!"]];
T[[All,k]]=DegradeM[(1/r[[k,k]])*T[[All,k]]];
];
{T,r}(*T=Q, r=R*)
]

Code for the Display of Results for errors in Decomposition

CompareDecomps[Good_,Bad_,A_,size_:450,log_:0,type_:1]:=Module[{Q,R,Q1,R1,labQR,l,n,m,tab={}},
{Q,R}=Good[A];
{Q1,R1}=Bad[A];
labQR={"Plot of Log|(Subscript[Q, Errors])-Q|","Plot of Log|(Subscript[R, Errors])-R|","Plot of Log|A-Q.R|","Plot of Log|A-Q.R|;Errors","Plot of Log|(Q^T)A-R|","Plot of Log|(Q^T)A-R|;Errors","Plot of I-Q.(Q^T); Correct Q","Plot of I-Q.(Q^T); With Errors","Histogram of Errors in A-Q.R","Histogram of Errors in A-Q.R;Degraded","Histogram of Errors in Log|A-Q.R|","Histogram of Errors in Log|A-Q.R|;Degraded","Histogram of Errors in (Q^T)A-R","Histogram of Errors in (Q^T)A-R;Degraded","Histogram of Errors in Log|(Q^T)A-R|","Histogram of Errors in Log|(Q^T)A-R|;Degraded"};
If[type==1,l=labQR;];
{m,n}=Dimensions[Q];

AppendTo[tab,
{MatrixPlot[Log[10,Abs[Q1-Q]],PlotLabel->l[[1]],ImageSize->size],
MatrixPlot[Log[10,Abs[R1-R]],PlotLabel->l[[2]],ImageSize->size]
}];

AppendTo[tab,
{MatrixPlot[Log[10,Abs[A-Q.R]],PlotLabel->l[[3]],ImageSize->size],
MatrixPlot[Log[10,Abs[A-Q1.R1]],PlotLabel->l[[4]],ImageSize->size]
}];

AppendTo[tab,
{MatrixPlot[Log[10,Abs[Transpose[Q].A-R]],PlotLabel->l[[5]],ImageSize->size],
MatrixPlot[Log[10,Abs[Transpose[Q1].A-R1]],PlotLabel->l[[6]],ImageSize->size]
}];

AppendTo[tab,
{MatrixPlot[Log[10,Abs[N[IdentityMatrix[m]-Q.Transpose[Q]]]],PlotLabel->l[[7]],ImageSize->size],
MatrixPlot[Log[10,Abs[N[IdentityMatrix[n]-Q1.Transpose[Q1]]]],PlotLabel->l[[8]],ImageSize->size]
}];
If[log==1||log==0,
AppendTo[tab,
{Histogram[Flatten[A-Q.R],PlotLabel->l[[9]],ImageSize->size],
Histogram[Flatten[A-Q1.R1],PlotLabel->l[[10]],ImageSize->size]
}];
];
If[log==2||log==0,
AppendTo[tab,
{Histogram[Log[10,Abs[Flatten[A-Q.R]]],PlotLabel->l[[11]],ImageSize->size],
Histogram[Log[10,Abs[Flatten[A-Q1.R1]]],PlotLabel->l[[12]],ImageSize->size]
}];
];
If[log==1||log==0,
AppendTo[tab,
{Histogram[Flatten[Transpose[Q].A-R],PlotLabel->l[[13]],ImageSize->size],
Histogram[Flatten[Transpose[Q1].A-R1],PlotLabel->l[[14]],ImageSize->size]
}];
];
If[log==2||log==0,
AppendTo[tab,
{Histogram[Log[10,Abs[Flatten[Transpose[Q].A-R]]],PlotLabel->l[[15]],ImageSize->size],
Histogram[Log[10,Abs[Flatten[Transpose[Q1].A-R1]]],PlotLabel->l[[16]],ImageSize->size]
}];
];

TableForm[tab,TableDirections->Column]
];

SolvR[R_,b_]:=Module[{y},
y=LinearSolve[Transpose[R],b];
LinearSolve[R,y]
];
SolvLU[L_,U_,b_]:=Module[{y},
y=LinearSolve[L,b];
LinearSolve[U,y]
];

Distributions Errors and Q’s (side topic)

Why do the Q’s look Normally distributed? There isn’t any apparent reason for this.

n=100;
A=GenMatrixConditionedAs[{n,n},10];
{Q,R}=QR[A];
{Q1,R1}=QRDecomposition[A];
K=EstimatedDistribution[Flatten[Q],NormalDistribution[mu,sig]];

size=600;
Show[
Histogram[Flatten[Q],Automatic,"PDF",ImageSize->size,PlotLabel->"Q's & Fitted Normal Dist"],
Plot[PDF[K,x],{x,-.5,.5},PlotStyle->Thick]
]
QuantilePlot[Flatten[Q],K,ImageSize->size,PlotLabel->"Q's vs Estimated Normal Dist"]
PearsonChiSquareTest[Flatten[Q]]

n=100;
reps=3000;
Timing[
pvalues=Table[
A=GenMatrixConditionedAs[{n,n},10];
{Q,R}=QR[A];
PearsonChiSquareTest[Flatten[Q]]
,{i,reps}];
];
Histogram[pvalues,Automatic,"PDF",ImageSize->size,PlotLabel->"Histogram of P-values of Chi-square tests"]
Print["P-value of Chi-Square Test checking p-values look Uniform: p = ",PearsonChiSquareTest[pvalues,UniformDistribution[{0,1}]]];

"P-value of Chi-Square Test checking p-values look Uniform: p = "0.00003188195768665381`

n=100;
reps=3000;
Timing[
pvalues=Table[
A=GenMatrixConditionedAs[{n,n},10];
{Q,R}=QR[A];
PearsonChiSquareTest[Flatten[Q]]
,{i,reps}];
];
Histogram[pvalues,Automatic,"PDF",ImageSize->size,PlotLabel->"Histogram of P-values of Chi-square tests"]
Print["P-value of Chi-Square Test checking p-values look Uniform: p = ",PearsonChiSquareTest[pvalues,UniformDistribution[{0,1}]]];

So, the distribution of the Q’s is non-normal because the p-values are unlikely to be uniform.

n=100;
reps=10000;
size=600;
Timing[ps=Table[
A=GenMatrixConditionedAsN[{n,n},10];
B=RandomReal[NormalDistribution[],{n,n}];
{Qa,Ra}=QR[A];
{Qb,Rb}=QR[B];
{PearsonChiSquareTest[Flatten[Qa]],PearsonChiSquareTest[Flatten[Qb]]}
,{i,reps}];]
Histogram[ps[[All,1]],Automatic,"PDF",ImageSize->size,PlotLabel->"Histogram of P-values of Chi-square tests"]
Print["P-value of Chi-Square Test checking p-values look Uniform: p = ",PearsonChiSquareTest[ps[[All,1]],UniformDistribution[{0,1}]]];
Histogram[ps[[All,2]],Automatic,"PDF",ImageSize->size,PlotLabel->"Histogram of P-values of Chi-square tests"]
Print["P-value of Chi-Square Test checking p-values look Uniform: p = ",PearsonChiSquareTest[ps[[All,2]],UniformDistribution[{0,1}]]];

n=300;
A=GenMatrixConditionedAs[{n,n},10];
{Q,R}=GSOn[A];
Histogram[Q[[All,k]]]
Histogram[Q[[k,All]]]
Manipulate[
Histogram[Flatten[Q[[All,k]]]]
,{k,1,n,1}];

Why don’t the distributions of the errors when the decomposition is put back together look normal? (They should) Those errors are just the result of a bunch of floating point computations, they should be normal.

n=100;
Zero1=A-Q.R;
Zero2=IdentityMatrix[n]-Q.Transpose[Q];
Norm[Zero1]
Norm[Zero2];

MatrixPlot[Log[10,Abs[Zero1]]];
MatrixPlot[Log[10,Abs[Zero2]]];

Histogram[-Log[10,Flatten[Zero1]]]
QuantilePlot[-Log[10,Flatten[Zero1]],ChiSquareDistribution[14.5]]
MatrixPlot[Log[Abs[Zero1]]];

MatrixPlot[Log[10,Abs[A-Q1.R1]]];

MatrixPlot[Log[10,Abs[A-Q.R]]];

Is there a reasonnig for the patterning of the errors when the decomposition is put back together?

ALSO: explore the same questions for the LU decomposition.

Propogation of Errors in Decompositions

ns=Table[
n=100;
A= GenMatrixConditionedAs[{n,n},n,1,1];
x=RandomReal[{-1,1},n];
b=A.x;
(*{L,U}=LU[A];*)
{Q,R}=GSO7[A];
x1=LinearSolve[R,Transpose[Q].b];
Norm[x-x1]/Norm[x]
,{20}];
Mean[ns]

What are the effects of the propagation of an error in the Cholesky Decomposition?

n=300;
A=GenMatrixConditionedAs[{n,n},n,1,1];
R=Chol[A];
R1=Cholloss[A,1];
R2=Cholloss[A,2];
Norm[Transpose[R].R-A]
Norm[Transpose[R1].R1-A]
Norm[Transpose[R2].R2-A]
MatrixPlot[R-R1];
MatrixPlot[R-R2];
MatrixPlot[Transpose[R1].R1-A];
MatrixPlot[Transpose[R2].R2-A];
x=RandomReal[{-1,1},n];
b=A.x;
xi=SolvR[R,b];
x1=SolvR[R1,b];
x2=SolvR[R2,b];
Print["Sol's norms"];
Norm[x-xi]
Norm[x-x1]
Norm[x-x2]
ListLogPlot[Abs[xi-x]];
ListLogPlot[Abs[x1-x]];
ListLogPlot[Abs[x2-x]];

n=100;
A=GenMatrixConditionedAs[{n,n},n,1,1];
{L,U}=LU[A];
{L1,U1}=LUloss1[A];
{L2,U2}=LUloss2[A];
(*
{L2,U2}=LUloss2[A];
{L3,U3}=LUloss3[A];
{L4,U4}=LUloss4[A];
{L5,U5}=LUloss5[A];
*)
Norm[L.U-A]
Norm[L1.U1-A]
Norm[L2.U2-A]
(*
Norm[L2.U2-A]
Norm[L3.U3-A]
Norm[L4.U4-A]
Norm[L5.U5-A]
*)
MatrixPlot[L1.U1-A]

What are the effects of the propagation of an error in the LU Decomposition?

What are the effects of the propagation of an error in the QR Decomposition?

A=GenMatrixConditionedAs[{100,100},100];
CompareDecomps[GSO,GSO5,A,250]

What are the effects of the propagation of an error in the GS Decomposition?

Work...

m=50;
ClearAll[A,Q,R,Q1,R1]
A=RandomReal[{-1,1},{m,m}];
Id=IdentityMatrix[m];
{Q,R}=GSO[A];
Norm[A-Q.R];
(*{Q2,R2}=GSO2[A];
Norm[A-Q2.R2]
Norm[A-(Q2^\[Transpose]).R2](*checking for orthogonality errors 8.07*)
Norm[Q-Q2]
Norm[R-R2]
Norm[Id-Q2.(Q2^\[Transpose])]
Norm[Id-(Q2^\[Transpose]).Q2]*)
(*____________________________________________________*)
{Q3,R3}=GSO3[A];
Norm[A-Q3.R3](*this gives the largest error*)
Norm[(Q3^\[Transpose]).A-R3](*checking for orthogonality errors*)
Norm[Q-Q3]
Norm[R-R3]
Norm[Id-Q3.(Q3^\[Transpose])]
(*{Q4,R4}=GSO4[A];
Norm[A-Q4.R4]
{Q5,R5}=GSO5[A];
Norm[A-Q5.R5]
{Q6,R6}=GSO6[A];
Norm[A-Q6.R6]*)
(*_________________________________________________
{Q7,R7}=GSO7[A];
Norm[A-Q7.R7](*followed by this one*)
Norm[(Q7^\[Transpose]).A-R7]
Norm[Q-Q7]
Norm[R-R7]
Norm[Id-Q7.(Q7^\[Transpose])](*checking for orthogonality errors*)*)

size=300;
TableForm[{
{MatrixPlot[Log[10,A-Q.R],PlotLabel->"Plot of Log|A-Q.R|",ImageSize->size],
MatrixPlot[Log[10,A-Q3.R3],PlotLabel->"Plot of Log|A-Q.R|;Errors",ImageSize->size]
},
{MatrixPlot[Log[10,(Q^\[Transpose]).A-R],PlotLabel->"Plot of Log|(Q^T)A-R|",ImageSize->size],
MatrixPlot[Log[10,(Q3^\[Transpose]).A-R3],PlotLabel->"Plot of Log|(Q^T)A-R|;Errors",ImageSize->size]
},
{MatrixPlot[Log[10,Q3-Q],PlotLabel->"Plot of Log|(Subscript[Q, Errors])-Q|",ImageSize->size],
MatrixPlot[Log[10,R3-R],PlotLabel->"Plot of Log|(Subscript[R, Errors])-R|",ImageSize->size]
},
{MatrixPlot[Log[10,Id-Q3.(Q3^\[Transpose])],PlotLabel->"Plot of I-Q.(Q^T); Correct Q",ImageSize->size],
MatrixPlot[Log[10,Id-Q3.(Q3^\[Transpose])],PlotLabel->"Plot of I-Q.(Q^T); With Errors",ImageSize->size]}
},TableDirections->Column]

Second Testing

ClearAll[G,S]
size=300;
{G,S}=GSOM[A];
TableForm[{
{MatrixPlot[Log[10,A-Q.R],PlotLabel->"Plot of Log|A-Q.R|",ImageSize->size],
MatrixPlot[Log[10,A-G.S],PlotLabel->"Plot of Log|A-Q.R|;Errors",ImageSize->size]
},
{MatrixPlot[Log[10,(Q^\[Transpose]).A-R],PlotLabel->"Plot of Log|(Q^T)A-R|",ImageSize->size],
MatrixPlot[Log[10,(G^\[Transpose]).A-S],PlotLabel->"Plot of Log|(Q^T)A-R|;Errors",ImageSize->size]
},
{MatrixPlot[Log[10,G-Q],PlotLabel->"Plot of Log|(Subscript[Q, Errors])-Q|",ImageSize->size],
MatrixPlot[Log[10,S-R],PlotLabel->"Plot of Log|(Subscript[R, Errors])-R|",ImageSize->size]
},
{MatrixPlot[Log[10,Id-Q.(Q^\[Transpose])],PlotLabel->"Plot of I-Q.(Q^T); Correct Q",ImageSize->size],
MatrixPlot[Log[10,Id-G.(G^\[Transpose])],PlotLabel->"Plot of I-Q.(Q^T); With Errors",ImageSize->size]}
},TableDirections->Column]

Change of error insertion point

Clear[data];
B=RandomReal[UniformDistribution[],{50,50}];
OrthB=Orthogonalize[B];(*Orthogonalizing B*)
data=Flatten[OrthB];
Histogram[Flatten[B]]
Histogram[data]
DistributionFitTest[Flatten[B],UniformDistribution[],Automatic]
ℋ=DistributionFitTest[data,UniformDistribution[],"HypothesisTestData"];
ℋ["TestDataTable",All]
(*QuantilePlot[data,ℋ["FittedDistribution"]]
DistributionFitTest[data,UniformDistribution[],Automatic]*)

Distribution of errors

ClearAll[A,Q,R,B,Q1,R1,data]
A=RandomReal[NormalDistribution[],{50,50}];
B=RandomReal[{0,1},{50,50}];
{Q,R}=GSO[A];
{Q1,R1}=GSO[B];
err=A-Q.R;
err1=B-Q1.R1;
data=Flatten[err];
data1=Flatten[err1];
Histogram[data]
Histogram[data1]
Norm[err]
Norm[err1]
Norm[IdentityMatrix[50]-Q.(Q^\[Transpose])]

Project Presentation #1

Code for Decompositions

Householder[x_List,i_]:= Module[{n=Length[x],u,H, τ=Norm[x[[i;;-1]]]},
H=SparseArray[Band[{1,1}]->1,{n,n}];
τ=Sign[x[[i]]] τ; 
u=Join[Table[0,{i-1}],{1},((x[[i+1;;-1]])/(τ+x[[i]]))];
H-(2/((Norm[u])^2))Outer[Times,u,u]
];
HouseholderDense[x_List,i_]:= Module[{n=Length[x],u,H, τ=Norm[x[[i;;-1]]]},
H=IdentityMatrix[n];
τ=Sign[x[[i]]] τ; 
u=Join[Table[0,{i-1}],{1},((x[[i+1;;-1]])/(τ+x[[i]]))];
H-(2/((Norm[u])^2))Outer[Times,u,u]
];

QR[A_List]:=Module[{R=A,n,m,Q,Qtemp,k},
{n,m}=Dimensions[A];
Q=IdentityMatrix[n];
If[n>m,k=m,k=n-1];
Do[
Qtemp=Householder[R[[All,j]],j];
R=Qtemp.R;
Q=Q.Qtemp;
,{j,1,k}];
{Q,R}
];

(*Gram-Schmidt with Double orthogonalization*)
GSO[A_]:=
Module[
{T=A,r,rHat,k,i,j,m=Length[(A^\[Transpose])]},
r=RandomReal[{-1,1},{m,m}]*0;(*Initializing r*)
For[k=1,k≤m,k++,
For[i=1,i≤k-1,i++,
r[[i,k]]=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*r[[i,k]](*first orthogonalization*);
];
For[i=1,i≤k-1,i++,
rHat=T[[All,k]].T[[All,i]];
T[[All,k]]-=T[[All,i]]*rHat(*second orthogonalization*);
r[[i,k]]+=rHat;
];
r[[k,k]]=Norm[T[[All,k]]];
If[r[[k,k]]≤$MachineEpsilon,
Return["The Matrix does not have a full rank!"]];
T[[All,k]]=(1/r[[k,k]])*T[[All,k]];
];
{T,r}(*T=Q, r=R*)
];

LU[A_]:=Module[{LU=A, n=Length[A],d,e,r,i,j},
Do[
d=1/LU[[r,r]];
Do[
e = d LU[[i,r]];
LU[[i,r]]=e;
Do[
LU[[i,j]]-= e* LU[[r,j]],
{j,r+1,n}
],
{i,r+1,n,1}],
{r,1,n-1,1}];
{LowerTriangularize[LU,-1]+SparseArray[Band[{1,1}]->1,{n,n}], 
UpperTriangularize[LU]}
];

The Propagation of Errors in the QR Decomposition

A=GenMatrixConditionedAs[{100,100},100];
CompareDecomps[QR,QRloss50x50,A,250]

A=GenMatrixConditionedAs[{100,100},100];
CompareDecomps[QR,QRloss1,A,250]

The Propagation of Errors in the GS Decomposition

A=GenMatrixConditionedAs[{100,100},100];
CompareDecomps[GSO,GSOM,A,300,2]

A=GenMatrixConditionedAs[{100,100},100];
CompareDecomps[GSO,GSO3,A,300,2]

The Propagation of Errors in the LU Decomposition

A=GenMatrixConditionedAs[{100,100},100];
CompareDecomps[LU,LUloss1,A,300]

A=GenMatrixConditionedAs[{60,60},100];
CompareDecomps[LU,LUloss2,A,300]

Distribution of the entries of Q

size=600;
Show[
Histogram[Flatten[Q],Automatic,"PDF",ImageSize->size,PlotLabel->"Q's & Fitted Normal Dist"],
Plot[PDF[K,x],{x,-.5,.5},PlotStyle->Thick]
]
QuantilePlot[Flatten[Q],K,ImageSize->size,PlotLabel->"Q's vs Estimated Normal Dist"]
PearsonChiSquareTest[Flatten[Q]]

n=100;
reps=3000;
Timing[
pvalues=Table[
A=GenMatrixConditionedAs[{n,n},10];
{Q,R}=QR[A];
PearsonChiSquareTest[Flatten[Q]]
,{i,reps}];
];
Histogram[pvalues,Automatic,"PDF",ImageSize->size,PlotLabel->"Histogram of P-values of Chi-square tests"]
Print["P-value of Chi-Square Test checking p-values look Uniform: p = ",PearsonChiSquareTest[pvalues,UniformDistribution[{0,1}]]];

"P-value of Chi-Square Test checking p-values look Uniform: p = "0.00003188195768665381`

n=100;
reps=3000;
Timing[
pvalues=Table[
A=GenMatrixConditionedAs[{n,n},10];
{Q,R}=QR[A];
PearsonChiSquareTest[Flatten[Q]]
,{i,reps}];
];
Histogram[pvalues,Automatic,"PDF",ImageSize->size,PlotLabel->"Histogram of P-values of Chi-square tests"]
Print["P-value of Chi-Square Test checking p-values look Uniform: p = ",PearsonChiSquareTest[pvalues,UniformDistribution[{0,1}]]];

So, the distribution of the Q’s is non-normal because the p-values are unlikely to be uniform.

n=100;
reps=10000;
size=600;
Timing[ps=Table[
A=GenMatrixConditionedAsN[{n,n},10];
B=RandomReal[NormalDistribution[],{n,n}];
{Qa,Ra}=QR[A];
{Qb,Rb}=QR[B];
{PearsonChiSquareTest[Flatten[Qa]],PearsonChiSquareTest[Flatten[Qb]]}
,{i,reps}];]
Histogram[ps[[All,1]],Automatic,"PDF",ImageSize->size,PlotLabel->"Histogram of P-values of Chi-square tests"]
Print["P-value of Chi-Square Test checking p-values look Uniform: p = ",PearsonChiSquareTest[ps[[All,1]],UniformDistribution[{0,1}]]];
Histogram[ps[[All,2]],Automatic,"PDF",ImageSize->size,PlotLabel->"Histogram of P-values of Chi-square tests"]
Print["P-value of Chi-Square Test checking p-values look Uniform: p = ",PearsonChiSquareTest[ps[[All,2]],UniformDistribution[{0,1}]]];

Some additional scratch work

Clear[B,data,Q,R,ℋ];
B=RandomReal[NormalDistribution[],{50,50}];
{Q,R}=GSO[B];
data=Flatten[Q];
(*Norm[IdentityMatrix[50]-Q.(Q^\[Transpose])]*)
Histogram[Flatten[B],PlotLabel->"Histogram showing the distribution of the entries in the Test matrix"]
Print["p=",DistributionFitTest[Flatten[B],Automatic],"\np>0.05, hence Normally distributed"]

Histogram[data,PlotLabel->"Histogram showing the distribution of the entries in Q"]

(*Test for goodness of fit*)

Print["p=",DistributionFitTest[data,Automatic],"\np>0.05, hence Normally distributed"]

ℋ=DistributionFitTest[data,Automatic,"HypothesisTestData"];
ℋ["TestDataTable",All]
Show[Histogram[data,Automatic,"ProbabilityDensity"],
Plot[PDF[ℋ["FittedDistribution"],x],{x,-1,1},PlotStyle->Thick]]
Print["There is insufficient evidence to reject a good fit to a Normal Distribution"]

DistributionFitTest[data,NormalDistribution[],{"TestDataTable","CramerVonMises"}]
QuantilePlot[data,ℋ["FittedDistribution"],PlotLabel->"Quantile plot- Verification of results"]

Timing[
j=Table[
n=300;
A=GenMatrixConditionedAs[{n,n},10];
{Q,R}=QR[A];
Correlation[Q[[All,n-1;;n]]][[2,1]]
,{20}];
Mean[j]
]

n=300;
A=GenMatrixConditionedAs[{n,n},10];
{Q,R}=QR[A];
Correlation[Q[[All,n-1;;n]]][[2,1]]

B

n=300;
A=GenMatrixConditionedAs[{n,n},n];
{L,U}=LU[A];
L[[60,50]]=Degrade2[L[[60,50]]];
x=RandomReal[{-1,1},n];
b=A.x;
x1=SolvLU[L,U,b];
ListPlot[x-x1]