Quick guide to Mathematica

To get help, place the cursor on the function and press F1
To run all commands in one cell, place the cursor there and press Shift + Enter or Enter on the numerical keyboard.
Output of many commands is suppressed by ; at the end of the command. If you want to see the output, delete ;
Argument of functions must be in square brackets [...], braces {...} are for arrays, ranges etc., use double brackets [[...]] to access elements of arrays
expr /. {x -> a} means substitute a for x in expr, expr // function means apply function to expr, i.e. it is equivalent to function[expr].
D[f[x],x] = derivative of f[x] with respect to x.
To define a function of x one can use f[x_] := 1+(x^2), notice the underscore and the colon, or # nad & like in f := 1+((#)^2) &
N[expr] evaluates expr with machine precision, if you use the decimal point in the expression it will be also evaluated with machine precision
Special characters can be inserted by pressing Esc ... Esc, e.g. Esc p Esc gives π
Clear[...] is used to unset any variables which could have been assigned previously.
Some other useful commands: Simplify, Expand, Factor; Integrate, Series, Sum

Clear all symbols from previous evaluations to avoid problems

Clear["Global`*"]

Examples of QR factorization

Run Basic.Functions.for.TMF057.m first ! To see more than 6 digits go to
Edit > Preferences > Appearance > Numbers > Formatting > Displayed precision 
and set 10 in Number of digits displayed in output.

Numerical loss of orthogonality using Gram-Schmidt orthogonalization

Testing matrix:

ϵ=(10^(-8));
tA={{1.0,1.0,1.0},{1.0+ϵ,1.0,1.0},{1.0+ϵ,1.0+ϵ,1.0}};
tA//MatrixForm

1) Standard Gram-Schmidt orthogonalization:

A=tA;
{Q,R}=MyQRGramSchmidt[A];
{Q//MatrixForm,R//MatrixForm}
Print["This should be an identity matrix:"];
Transpose[Conjugate[Q]].Q//MatrixForm

2) Modified Gram-Schmidt orthogonalization:

A=tA;
{Q,R}=MyQRModifiedGramSchmidt[A];
{Q//MatrixForm,R//MatrixForm}
Print["This should be an identity matrix:"];
Transpose[Conjugate[Q]].Q//MatrixForm

3) Householder triangularization, returns (SuperDagger[Q]):

A=tA;
{Q,R}=MyQRHouseholder[A];
{Transpose[Q]//MatrixForm,R//MatrixForm}
Print["This should be an identity matrix:"];
Transpose[Conjugate[Q]].Q//MatrixForm

4) Mathematica built-in algorithm - also Householder triangularization, but probably with column pivoting (even more precise):

A=tA;
{Q,R}=QRDecomposition[A];
{Transpose[Q]//MatrixForm,R//MatrixForm}
Print["This should be an identity matrix:"];
Transpose[Conjugate[Q]].Q//MatrixForm

Instability of basic Gram-Schmidt orthogonalization

m=101;
x = Range[-1.0,1.0,2.0/(m-1)];
A=ConstantArray[1.0,{m,m}];
A[[All,2;;m]]=Table[Table[x[[im]]^in,{in,1,m-1}],{im,1,m}]; (* Vandermode matrix *)
{Q1,R1}=MyQRGramSchmidt[A];
{Q2,R2}=MyQRModifiedGramSchmidt[A];
{Q3,R3}=MyQRHouseholder[A];
{Q4,R4}=QRDecomposition[A];

ListLogPlot[{Diagonal[R1],Diagonal[R2],Diagonal[R3],Diagonal[R4]},PlotStyle->{Red,Blue,Green,Orange}]

Backward stability of Householder triangularization

Prepare a matrix A as Q.R where Q and R are random matrices:

n=50;
R=UpperTriangularize[RandomReal[{0,1},{n,n}]];
{Q,Rtemp}=QRDecomposition[RandomReal[{-1,1},{n,n}]];
A=Q.R;

Calculate QR factorization. How accurate are calculated Q1 and R1?

{Q1,R1}=MyQRHouseholder[A];
Print["Error in Q: ",Norm[ConjugateTranspose[Q1]-Q]/Norm[Q]];
Print["Error in R: ",Norm[R1-R]/Norm[R]];
{U,S,V}=SingularValueDecomposition[A];
σ=Diagonal[S];
κA=Max[σ]/Min[σ];
Print["Condition number: ",κA]

How accurate is Q1.R1?

Norm[A-ConjugateTranspose[Q1].R1]/Norm[A]