(1.1) Normal equation: A^*Ax=A^*b,

where, matrix A^*A:
184.9198 182.5428 165.9724 122.5413
182.5428 273.5335 199.8720 157.2537
165.9724 199.8720 268.2798 216.0125
122.5413 157.2537 216.0125 232.1604

Matrix A^*b:
98.0007
175.8608
182.8303
193.3366

(1.2) C-source code of Cholesky decomposition method:

#include "stdafx.h"
#include "math.h"


int _tmain(int argc, _TCHAR* argv[])
{
int i, j, k, m;

double aa[8][4]={7.0605, 9.5022, 4.4559, 1.19,
0.3183, 0.3445, 6.4631, 4.9836,
2.7692, 4.3874, 7.0936, 9.5974,
0.4617, 3.8156, 7.5469, 3.4039,
0.9713, 7.6552, 2.7603, 5.8527,
8.2346, 7.952, 6.797, 2.2381,
6.9483, 1.8687, 6.551, 7.5127,
3.171, 4.8976, 1.6261, 2.551};

double bb[8] = {0.7585,
0.5395,
5.3080,
7.7917,
9.3401,
1.2991,
5.6882,
4.6939};

double l[4][4];
double y[4];
double x[4];
double AC[4][4];
double ax[4];
double xx[4];
double at[4][8];
double a[4][4];
double b[4];


for(i = 0; i<4; i ++){
for(j = 0; j<8 ; j ++){
at[i][j]=0;
at[i][j]=aa[j][i];
}
}

for(i=0; i<4; i++){
for(j=0; j<4; j++){
a[i][j]=0;
for (k = 0; k<8; k++) {
a[i][j]+=at[i][k]*aa[k][j];
}
}
}
printf("Matrix At*A: \n");
for(i=0; i<4; i++){
for(j=0; j<4; j++){
printf(" %7.4f", a[i][j]);
}
printf("\n");
}

for(i=0; i<4; i++){
b[i]=0;
for (k = 0; k < 8; k++) {
b[i]+=at[i][k]*bb[k];
}
}
printf("Matrix At*b: \n");
for(i=0; i<4; i++){
printf(" %7.4f", b[i]);
printf("\n");
}


for (i = 0; i <4; i++){
for (j = 0; j <4; j++){
l[i][j]=0;
}
}
for (i = 0; i <4; i++){
y[i]=b[i];
}
for(i=0; i<4; i++){
for(j=0; j<4; j++){
AC[i][j]=0;
AC[i][j]=a[i][j];
}
}

/* Cholesky Decomposition */
for (k = 0; k < 4-1; k++){
l[k][k]=sqrt(a[k][k]);
for(m=k+1; m<4; m++){
l[m][k]=a[m][k]/l[k][k];
for(j=k+1; j<4; j++){
for(i=m; i<4; i++){
a[i][j]-=l[i][k]*l[j][k];
}
}
}
}
l[3][3]=sqrt(a[3][3]);

printf("Matrix B: \n");
for(i=0; i<4; i++){
for(j=0; j<4; j++){
printf(" %7.4f", l[i][j]);
}
printf("\n");
}

for (i = 0; i <4; i++) {
y[i] /= l[i][i];
for (j = 0; j <i; j++){
y[i] -= l[i][j] * y[j]/l[i][i];
}
}
for (i = 0; i <4; i++){
x[i]=y[i];
}
for (i =4-1; i >=0; i--) {
x[i] /= l[i][i];
for (j = i+1; j < 4; j++){
x[i] -= l[j][i] * x[j]/l[i][i];
}
}

printf("x= \n");
for(i=0; i<4; i++){
printf(" %7.4f", x[i]);
printf("\n");
}

for(i=0; i<4; i++){
ax[i]=0;
for (k = 0; k < 4; k++) {
ax[i]+=AC[i][k]*x[k];
}
}
for(i=0; i<4; i++){
xx[i]=ax[i]-b[i];
}
printf("At*Ax-At*b= \n");
for(i=0; i<4; i++){
printf(" %7.4f", xx[i]);
printf("\n");
}

return 0;
}
 

Using Cholesky decomposition method, we obtain A=B^tB,

where, matrix B:
13.5985 0.0000 0.0000 0.0000
13.4237 9.6611 0.0000 0.0000
12.2052 3.7297 10.2666 0.0000
9.0114 3.7561 8.9629 7.5176
 

and the solution is

x=
-0.4958
0.5741
-0.0304
0.7339
 

We check the solution is right using the residual:

At*Ax-At*b=
0.0000
0.0000
-0.0000
0.0000

(2.1) Find the Householder QR factorization for

First of all, we find the Householder matrix for the transformation of the first column vector x = of A.

we compute Choose the sign plus for coefficient of , and we obtain    Find the Householder matrix H₁ that depends only on direction of v,

Thus,

We can verify that

To find we compute the according Householder vector: the first column vector of the (1,1)-minor of H₁A:

Hence

and

and also



Find also the orthogonal matrix



and check the result

(2.2) Find the Householder QR factorization of


The vector that has to be transformed is where Construct the vector

Choose the sign minus for the coefficient of and take into account that H depends only on the direction of :

Find the Householder matrix



Verify that H₁ annihilates all the elements of the first column of A but the first one. Indeed,

Further we transform the vector the first column vector of the (1,1)-minor of H₁A, where

Find the Householder vector according to

Choose the sign minus for coefficient of :

We obtain the Householder matrix according to this vector





and find that

Thus,

and

Let us check the result: