Linear Algebra and the C Language/a0l9

__NOTOC__ 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.