Linear Algebra and the C Language/a0l9
Install and compile this file in your working directory.
/* ------------------------------------ */
/* Save as: c00c2.c */
/* ------------------------------------ */
#include "v_a.h"
/* ------------------------------------ */
/* ------------------------------------ */
#define RA R3
#define CA C2
/* ------------------------------------ */
/* ------------------------------------ */
int main(void)
{
double tA[RA*CA]={
+0.078326044999, -0.989057126856,
+0.704934404989, +0.143625457365,
+0.704934404989, -0.033730221048
};
double tx[RA*C1]={
+5.0000,
+8.0000,
-8.0000,
};
double **A = ca_A_mR(tA, i_mR(RA,CA));
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 = ca_A_mR(tx, i_mR(RA,C1));
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);
stop();
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);
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 R3
Find a transformation matrix for
a projection onto R3 :
Proj(x) = A inv(AT A) AT x
A:
+0.0783 -0.9891
+0.7049 +0.1436
+0.7049 -0.0337
Compute Proj(x) with:
x:
+5.0000
+8.0000
-8.0000
Press return to continue.
AT:
+0.0783 +0.7049 +0.7049
-0.9891 +0.1436 -0.0337
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.0783 +0.7049 +0.7049
-0.9891 +0.1436 -0.0337
V = A inv(AT A) AT:
+0.9844 -0.0868 +0.0886
-0.0868 +0.5176 +0.4921
+0.0886 +0.4921 +0.4981
Press return to continue.
V is transformation matrix for
a projection onto a R3 subspace
V:
+0.9844 -0.0868 +0.0886
-0.0868 +0.5176 +0.4921
+0.0886 +0.4921 +0.4981
inv( ID )
Proj(x) = A inv(AT A) AT x
Proj(x) = V x:
+3.5185
-0.2304
+0.3950
Press return to continue.
inv( ID )
Proj(x) = A inv(AT*A) AT x
+3.5185
-0.2304
+0.3950
A AT:
+0.9844 -0.0868 +0.0886
-0.0868 +0.5176 +0.4921
+0.0886 +0.4921 +0.4981
Proj(x) = (A AT) x
+3.5185
-0.2304
+0.3950
Press return to continue.