Linear Algebra and the C Language/a0l8
Install and compile this file in your working directory.
/* ------------------------------------ */
/* Save as: c00c1.c */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
/* ------------------------------------ */
#define RA R4
#define CA C2
/* ------------------------------------ */
/* ------------------------------------ */
void fun(void)
{
double **A = r_Q_mR( i_mR(RA,CA),9);
double **AT = transpose_mR(A, i_mR(CA,RA));
double **ATA = mul_mR(AT,A, i_mR(CA,CA));
double **invATA = invgj_mR(ATA, i_mR(CA,CA));
double **invATA_AT = mul_mR(invATA,AT, i_mR(CA,RA));
double **V = mul_mR(A,invATA_AT, i_mR(RA,RA));
double **AAT = mul_mR(A,AT, i_mR(RA,RA));
double **x = r_mR( i_mR(RA,C1),9.);
double **Vx = mul_mR(AAT,x, i_mR(RA,C1));
clrscrn();
printf(" A is subspace of R%d \n\n"
" Find a transformation matrix for \n"
" a projection onto R%d : \n\n"
" Proj(x) = A inv(AT A) AT x \n\n",RA,RA);
printf(" A:");
p_mR(A, S5,P4,C7);
printf(" Compute Proj(x) with: \n\n"
" x:");
p_mR(x, S5,P4,C7);
stop();
clrscrn();
printf(" AT:");
p_mR(AT, S5,P4,C7);
printf(" ATA:");
p_mR(ATA, S5,P4,C7);
printf(" inv(AT A):");
p_mR(invATA, S5,P4,C7);
printf(" inv(AT A) AT:");
p_mR(invATA_AT, S5,P4,C7);
printf(" V = A inv(AT A) AT:");
p_mR(V, S5,P4,C7);
stop();
clrscrn();
printf(" V is transformation matrix for \n"
" a projection onto a R%d subspace\n\n"
" V:",RA);
p_mR(V, S5,P4,C7);
printf(" inv( ID ) \n"
" Proj(x) = A inv(AT A) AT x\n\n");
printf(" Proj(x) = V x:");
p_mR(mul_mR(V,x,Vx), S5,P4,C7);
stop();
clrscrn();
printf(" inv( ID ) \n"
" Proj(x) = A inv(AT*A) AT x");
p_mR(Vx, S5,P4,C7);
printf(" A AT:");
p_mR(AAT, S5,P4,C7);
printf(" Proj(x) = (A AT) x");
p_mR(Vx, S5,P4,C7);
f_mR(A);
f_mR(AT);
f_mR(ATA);
f_mR(invATA);
f_mR(invATA_AT);
f_mR(V); /* A inv(AT A) AT */
f_mR(AAT);
f_mR(x);
f_mR(Vx);
}
/* ------------------------------------ */
int main(void)
{
time_t t;
srand(time(&t));
do
{
fun();
} while(stop_w());
return 0;
}
/* ------------------------------------ */
/* ------------------------------------ */
The orthonormality condition can also be dropped. If u_n is a (not necessarily orthonormal) basis with , and A is the matrix with these vectors as columns, then the projection is:
Proj(x) = A inv(AT A) AT
If A is orthonomal then (AT A) = ID and then:
Proj(x) = A inv(ID) AT= (A AT)
... Wikipedia: Projection (linear algebra)
Screen output example:
A is subspace of R4
Find a transformation matrix for
a projection onto R4 :
Proj(x) = A inv(AT A) AT x
A:
+0.7493 +0.1303
-0.4683 +0.5568
-0.2810 -0.8032
+0.3746 -0.1669
Compute Proj(x) with:
x:
+4.0000
+6.0000
-1.0000
+4.0000
Press return to continue.
AT:
+0.7493 -0.4683 -0.2810 +0.3746
+0.1303 +0.5568 -0.8032 -0.1669
ATA:
+1.0000 +0.0000
+0.0000 +1.0000
inv(AT A):
+1.0000 -0.0000
-0.0000 +1.0000
inv(AT A) AT:
+0.7493 -0.4683 -0.2810 +0.3746
+0.1303 +0.5568 -0.8032 -0.1669
V = A inv(AT A) AT:
+0.5784 -0.2783 -0.3152 +0.2589
-0.2783 +0.5294 -0.3157 -0.2684
-0.3152 -0.3157 +0.7240 +0.0288
+0.2589 -0.2684 +0.0288 +0.1682
Press return to continue.
V is transformation matrix for
a projection onto a R4 subspace
V:
+0.5784 -0.2783 -0.3152 +0.2589
-0.2783 +0.5294 -0.3157 -0.2684
-0.3152 -0.3157 +0.7240 +0.0288
+0.2589 -0.2684 +0.0288 +0.1682
inv( ID )
Proj(x) = A inv(AT A) AT x
Proj(x) = V x:
+1.9946
+1.3049
-3.7634
+0.0695
Press return to continue.
inv( ID )
Proj(x) = A inv(AT*A) AT x
+1.9946
+1.3049
-3.7634
+0.0695
A AT:
+0.5784 -0.2783 -0.3152 +0.2589
-0.2783 +0.5294 -0.3157 -0.2684
-0.3152 -0.3157 +0.7240 +0.0288
+0.2589 -0.2684 +0.0288 +0.1682
Proj(x) = (A AT) x
+1.9946
+1.3049
-3.7634
+0.0695
Press return to continue
Press X return to stop