Linear Algebra and the C Language/a0l8

__NOTOC__ 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