ill-posed Matrix and Condition Number

Example of an ill-posed system: Hilbert matrix

We consider the linear system of equations (Subscript[H, n ])x = b, in which the matrix (Subscript[H, n ]) is the Hilbert matrix, (Subscript[H, n ])= (Subscript[((1/(i+j+1))), i,j=0]^n), and we choose b at the right hand side so that the solution is given by x = (1,1,...,1).

H[n_]:=Table[1/(i+j+1), {i,0,n},{j,0,n}];
RHS[n_]:=Table[Sum[ 1/(i+1), {i,k,k+n}],{k,0,n}];
X0[n_]:= Table[1, {i,0,n}];

StringForm["H[6]= `` RHS[6]=``", MatrixForm[H[6]],  MatrixForm[RHS[6]]]

We solve this system of equations with the Mathematica solver NSolve. The solution is a vector X.

Var[n_]:=Table[x[i],{i,1,n+1}]
sol[n_]:=Solve[N[H[n]].Var[n]==N[RHS[n]]];
X[n_]:=Table[sol[n][[1,i,2]],{i,1,n+1}]

Here is the solution for the system with n = 8, 9,10, 11, 12, ...

Do[Print[StringForm["n=``,     X[n]= `` " , n , X[n] ]],   {n, 6, 12}]

Below is the residue vector for the case n =12

Max[Abs[H[12].X[12] - RHS[12] ]]

Below is the error vector for the case n =12

Max[Abs[X[12]-X0[12]]]

This example shows that the residue vector may not be a good measure for the accuracy of the solution. Indeed, in the above example, the residue vector is as small as 10^{-16}, but the solution is wrong (the error is > 5). The reason for the problem lies in the ill-condition of the Hilbert matrix. Notice that the Mathematica works with the machine accuracy = 16:

$MachinePrecision

Example of an ill-posed system: Vandermond matrix

We consider the linear system of equations (Subscript[V, n ])x = b, in which the matrix (Subscript[V, n ]) is the Vandermond matrix, (Subscript[V, n ])= (Subscript[((((i+1))^j)), i,j=0]^n), and we choose b at the right hand side so that the solution is given by x = (1,1,...,1).

VV[n_]:=Table[i^(j-1),{i,1,n},{j,1,n}];
RHS[n_]:=Table[Sum[VV[n][[j,i]],{i,1,n}],{j,1,n}];
X0[n_]:= Table[1, {i,1,n}];

StringForm["VV[n]= ``      RHS[n]=``.",MatrixForm[VV[6]],MatrixForm[RHS[6]]]

We solve this system of equations with the Mathematica solver NSolve. The solution is a vector X.

Var[n_]:=Table[x[i],{i,1,n}]
sol[n_]:=NSolve[N[VV[n]].Var[n]==N[RHS[n]]];
X[n_]:=Table[sol[n][[1,i,2]],{i,1,n}]

Here is the solution for the system with n = 13, 14,15, 16, 17, ...

Do[Print[StringForm["n=``,     X[n]= `` " , n , X[n] ]],   {n, 10, 15}]

Below is the residue vector for the case n =14

Max[Abs[VV[15].X[15] - RHS[15] ]]

Below is the error vector for the case n =14

Max[Abs[X[15]-X0[15]]]

This example shows that the residue vector may not be a good measure for the accuracy of the solution. Indeed, in the above example, the residue vector is as small as 10^{-16}, but the solution is wrong (the error is > 5). The reason for the problem lies in the ill-condition of the Vandermond matrix. Notice that the Mathematica works with the machine accuracy = 16:

$MachinePrecision

Condition numbers

Below is a program that computes the conditional number Cond(A) of a matrix A = (Subscript[((Subscript[a, i,j]) ), i,j=0]^n). The condition number is defined by
                                 Cond(A) = || A || × ||(A^(-1))||,          where     ||A|| = (max_(1≤i≤n))  (∑_(j=1)^n)|(Subscript[a, i,j])|.

Cond[H0_] := Module[{H, invH, n, h, g,x},
     H=H0; n=Length[H]; invH=Inverse[H]; Module[{i,j},h=Table[0,{n}]; g=Table[0,{n}];
     For[i=1,i<=n,i++, h[[i]]=Sum[Abs[H[[i,j]]],  {j,1,n}] ];  
     For[i=1,i<=n,i++, g[[i]]=Sum[Abs[invH[[i,j]]],  {j,1,n}] ]; ]; 
x= Max[Table[h[[i]], {i,1,n} ]]* Max[Table[g[[i]], {i,1,n} ]];  
Return[x]];

It is known that the condition number gauges the transfer of errors from the matrix A and the vector b to the solution of Ax = b. When the condition number is large, the matrix is called ill-conditioned and the system of equations is ill-posed.  The Hilbert matrix (Subscript[H, n]) is known to be severely ill-conditioned even for moderate n. Below we compute its condition number:

Do[Print[StringForm["n=``,     Cond[Hn]= `` " , n , Cond[H[n]]//N  ]],   {n, 8, 12}]

The correlation between the condition numbers and the errors of the solutions in the previous example is clearly visible.

H[n_]:=Table[1/(i+j+1), {i,0,n},{j,0,n}];
RHS[n_]:=Table[Sum[ 1/(i+1), {i,k,k+n}],{k,0,n}];
X0[n_]:= Table[1, {i,0,n}]; 
Var[n_]:=Table[x[i],{i,1,n+1}]
sol[n_]:=Solve[N[H[n]].Var[n]==N[RHS[n]]];
X[n_]:=Table[sol[n][[1,i,2]],{i,1,n+1}]

Do[Print[StringForm["n=``,     X[n]= `` " , n , X[n] ]],   {n, 8, 12}]