{
"cells": [
 {
		   "cell_type": "markdown",
	   "metadata": {},
	   "source": [
       "# Chapter 9: Linear Least Squares Problem"
	   ]
	},
{
		   "cell_type": "markdown",
		   "metadata": {},
		   "source": [
			"## Example 9.1: Moore_Penrose_Generalized_Inverse.sce"
		   ]
		  },
  {
"cell_type": "code",
	   "execution_count": null,
	   "metadata": {
	    "collapsed": true
	   },
	   "outputs": [],
"source": [
"//Example 9.1\n",
"//Moore-Penrose Generalized Inverse\n",
"//Page no. 292\n",
"clc;clear;close;\n",
"\n",
"AT=[3,0,3;0,3,3];\n",
"A=AT';           //transpose\n",
"I=inv(AT*A);       //inverse\n",
"disp(I,'Inverse of AT*A=',AT*A,'AT*A=',A,'A=',AT,'AT=');\n",
"A#=I*AT;\n",
"disp(A#,'Moore-Penrose Generalized Inverse of A=')"
   ]
   }
,
{
		   "cell_type": "markdown",
		   "metadata": {},
		   "source": [
			"## Example 9.2: Curve_Fitting.sce"
		   ]
		  },
  {
"cell_type": "code",
	   "execution_count": null,
	   "metadata": {
	    "collapsed": true
	   },
	   "outputs": [],
"source": [
"//Example 9.2\n",
"//Curve Fitting\n",
"//Page no. 293\n",
"clc;clear;close;\n",
"x(1)=0.25;\n",
"for i=2:6\n",
"    x(1,i)=x(1,i-1)+0.25;\n",
"end               //x values\n",
"y(1,1)=3.1;y(1,2)=1.7;y(1,3)=1;y(1,4)=0.68;y(1,5)=0.42;y(1,6)=0.26; //y values\n",
"\n",
"//construction of normal equations\n",
"for i=1:6\n",
"    Y(1,i)=log10(y(1,i));\n",
"end\n",
"Ex=0;\n",
"for i=1:6\n",
"    Ex=Ex+x(1,i);\n",
"end\n",
"EY=0;\n",
"for i=1:6\n",
"    EY=EY+Y(1,i);\n",
"end\n",
"Ex2=0;\n",
"for i=1:6\n",
"    Ex2=Ex2+x(1,i)^2;\n",
"end\n",
"ExY=0;\n",
"for i=1:6\n",
"    ExY=ExY+x(1,i)*Y(1,i);\n",
"end\n",
"printf('E x(k)\t\t y(k)\t\tE Y(k)\t\tE x2(k)\t\tE x(k)*Y(k)')\n",
"printf('\n----------------------------------------------------------------------------')\n",
"for i=1:6\n",
"    printf('\n%f\t%f\t%f\t%f\t%f',x(1,i),y(1,i),Y(1,i),x(1,i)^2,x(1,i)*Y(1,i))\n",
"end\n",
"printf('\n----------------------------------------------------------------------------')\n",
"printf('\n%f\t%f\t%f\t%f\t%f',Ex,0,EY,Ex2,ExY)\n",
"printf('\n----------------------------------------------------------------------------\n\n')\n",
"A=[6,Ex;Ex,Ex2];      //system of normal equations\n",
"B=[EY;ExY];\n",
"X=inv(A)*B;\n",
"a=exp(X(1,1));\n",
"b=-1*X(2,1);\n",
"for i=1:2\n",
"    for j=1:2\n",
"        printf('%f    ',A(i,j))\n",
"    end\n",
"    if(i==1)\n",
"        printf('  *')\n",
"    end\n",
"    \n",
"    printf('\ta%i',i);\n",
"    if(i==1)\n",
"        printf('  =')\n",
"    end\n",
"    \n",
"    printf('\t%f\n',B(i,1))\n",
"end\n",
"printf('\n\na1=%f\na2=%f\n\na=%f\nb=%f\n\n',X(1,1),X(2,1),a,b)\n",
"printf('The fitted curve is:\n            %fx\ny=%f e',b,a)"
   ]
   }
,
{
		   "cell_type": "markdown",
		   "metadata": {},
		   "source": [
			"## Example 9.3: Gram_Schmidt_Orthogonalization_or_Orthonormalization_Process.sce"
		   ]
		  },
  {
"cell_type": "code",
	   "execution_count": null,
	   "metadata": {
	    "collapsed": true
	   },
	   "outputs": [],
"source": [
"//Example 9.3\n",
"//Gram-Schmidt Orthogonalization/Orthonormalization Process\n",
"//Page no. 294\n",
"clc;clear;close;\n",
"deff('y=f(x,a)','y=sqrt(x(1,a)^2+x(2,a)^2+x(3,a)^2+x(4,a)^2)');\n",
"deff('y=f1(g,a,h,b)','y=g(1,a)*h(1,b)+g(2,a)*h(2,b)+g(3,a)*h(3,b)+g(4,a)*h(4,b)');\n",
"\n",
"U=[1/sqrt(3),-2/sqrt(7),1,0,0,0;0,1/sqrt(7),0,1,0,0;1/sqrt(3),1/sqrt(7),0,0,1,0;-1/sqrt(3),-1/sqrt(7),0,0,0,1];\n",
"for i=1:4\n",
"        V(i,1)=U(i,1);\n",
"end\n",
"for i=1:4\n",
"    if(f(V,1)~=0)\n",
"        W(i,1)=V(i,1)/f(V,1);\n",
"    else\n",
"        W(i,1)=0;\n",
"    end \n",
"end\n",
"for j=2:6\n",
"    for i=1:4\n",
"        for l=1:4\n",
"            k(l,1)=0;\n",
"        end\n",
"        for l=1:j-1\n",
"            for m=1:4\n",
"                w(m,1)=W(m,l);\n",
"            end\n",
"            k=k-(f1(U,j,W,l))*w;\n",
"        end\n",
"        V(i,j)=U(i,j)+k(i,1);\n",
"    end\n",
"    for i=1:4\n",
"        if(j~=4)\n",
"            if(f(V,j)~=0)\n",
"                W(i,j)=V(i,j)/f(V,j);\n",
"            else\n",
"                W(i,j)=0;\n",
"            end \n",
"        else\n",
"            W(i,j)=0;\n",
"        end\n",
"    end\n",
"    \n",
"end\n",
"disp(U,'U=')\n",
"disp('W=')\n",
"printf('\n')\n",
"for i=1:4\n",
"    for j=1:6\n",
"        printf('%.4f\t\t',W(i,j))\n",
"    end\n",
"    printf('\n')\n",
"end\n",
"disp('V=')\n",
"printf('\n')\n",
"for i=1:4\n",
"    for j=1:6\n",
"        printf('%.4f\t\t',V(i,j))\n",
"    end\n",
"    printf('\n')\n",
"end"
   ]
   }
,
{
		   "cell_type": "markdown",
		   "metadata": {},
		   "source": [
			"## Example 9.4: QR_Decompositio.sce"
		   ]
		  },
  {
"cell_type": "code",
	   "execution_count": null,
	   "metadata": {
	    "collapsed": true
	   },
	   "outputs": [],
"source": [
"//Example 9.4\n",
"//QR Decomposition\n",
"//Page no. 296\n",
"clc;clear;close;\n",
"\n",
"A=[2,1,1;1,3,1;1,1,4];\n",
"B=A*A';\n",
"disp(B,'AT*A=')\n",
"//cholesky factorization to find R\n",
"R(2,1)=0;R(3,1)=0;R(3,2)=0;     \n",
"for i=1:3\n",
"    for j=1:3\n",
"        if(i==j)\n",
"            k=0;\n",
"            for m=1:i-1\n",
"               k=k+R(m,i)^2; \n",
"            end\n",
"            R(i,j)=sqrt(B(i,j)-k)\n",
"        end\n",
"        if(j>i)\n",
"            k=0;\n",
"            for m=1:i-1\n",
"                k=k+R(m,j)*R(m,i);\n",
"            end\n",
"            R(i,j)=(B(i,j)-k)/R(i,i)\n",
"        end\n",
"    end\n",
"end\n",
"//cholesky factorization end\n",
"disp(R,'Upper Triangular Matrix (R)=')\n",
"R_1=inv(R);\n",
"disp(R_1,'Inverse of R')\n",
"Q=A*R_1;\n",
"disp(Q,'Orthogonal Matrix Q=')"
   ]
   }
,
{
		   "cell_type": "markdown",
		   "metadata": {},
		   "source": [
			"## Example 9.5: Vector_Computatio.sce"
		   ]
		  },
  {
"cell_type": "code",
	   "execution_count": null,
	   "metadata": {
	    "collapsed": true
	   },
	   "outputs": [],
"source": [
"//Example 9.5\n",
"//Vector Computation\n",
"//Page no. 299\n",
"clc;clear;close;\n",
"\n",
"X=[2,3,0,1];\n",
"n=X(1);\n",
"for i=2:4\n",
"    if(n<X(i))\n",
"        n=X(i);\n",
"    end\n",
"end\n",
"printf('\nMaximum Value (n)=%i\n',n)\n",
"for i=1:4\n",
"    X(i)=X(i)/n;\n",
"end\n",
"disp(X,'Normalized X=')\n",
"k=0;\n",
"for i=1:4\n",
"    k=k+X(i)^2;\n",
"end\n",
"sigma=X(1)*abs(1/X(1))*sqrt(k);\n",
"printf('\nsigma=%f\n',sigma);\n",
"X(1)=X(1)+sigma;\n",
"printf('\nModified x1 = %g\n',X(1))\n",
"for i=1:4\n",
"    U(1,i)=X(i);\n",
"end\n",
"disp(U,'U=')\n",
"p=sigma*X(1);sigma=n*sigma;\n",
"printf('\n p = %f\n\n sigma = %f',p,sigma);\n",
"printf('\n\n\nNote : There is a computation error in calculation of U1')"
   ]
   }
,
{
		   "cell_type": "markdown",
		   "metadata": {},
		   "source": [
			"## Example 9.6: House_Holder_Transformation.sce"
		   ]
		  },
  {
"cell_type": "code",
	   "execution_count": null,
	   "metadata": {
	    "collapsed": true
	   },
	   "outputs": [],
"source": [
"//Example 9.6\n",
"//House Holder Transformation\n",
"//Page no. 300\n",
"clc;clear;close;\n",
"\n",
"A=[4,2,1;2,5,-2;1,-2,7]\n",
"disp(A,'A=')\n",
"k=0;\n",
"for j=2:3\n",
"    k=k+A(j,1)^2;\n",
"end\n",
"a=A(2,1)*abs(1/A(2,1))*sqrt(k);\n",
"disp(a,'alpha=')\n",
"U=[0;a+A(2,1);A(3,1)];\n",
"disp(U,'U=')\n",
"U1=U'*U;\n",
"disp(U1,'UT*U=')\n",
"U2=U*U';\n",
"disp(U2,'U*UT=')\n",
"P=eye(3,3)-(2*U2)/U1;\n",
"disp(P,'P=');\n",
"B=P*A*P;\n",
"disp(B,'B=');\n",
"printf('\n\n\nThere are computation error in the answers given by the book in this example\n\n(a22 value error in U*UT)')"
   ]
   }
,
{
		   "cell_type": "markdown",
		   "metadata": {},
		   "source": [
			"## Example 9.7: Givens_QR_Method.sce"
		   ]
		  },
  {
"cell_type": "code",
	   "execution_count": null,
	   "metadata": {
	    "collapsed": true
	   },
	   "outputs": [],
"source": [
"//Example 9.7\n",
"//Givens QR Method\n",
"//Page no. 303\n",
"clc;clear;close;\n",
"\n",
"A=[4,2,1;2,5,-2;1,-2,7]\n",
"deff('y=c(i,j)','y=A(j,j)/sqrt((A(i,j)^2+A(j,j)^2))')\n",
"deff('y=s(i,j)','y=A(i,j)/sqrt((A(i,j)^2+A(j,j)^2))')\n",
"disp(A,'A=')\n",
"R=A;Q=eye(3,3);\n",
"m=1;\n",
"for j=1:2\n",
"    for i=j+1:3\n",
"        for k=1:3             //C matrix evaluation\n",
"            for l=1:3\n",
"                if(k==l)\n",
"                    if(k==i | k==j)\n",
"                        C(k,l)=c(i,j)\n",
"                    else\n",
"                        C(k,l)=1\n",
"                    end\n",
"                end\n",
"                if(k>l)\n",
"                    if(k==i & l==j)\n",
"                        C(k,l)=-1*s(i,j)\n",
"                    else\n",
"                        C(k,l)=0\n",
"                    end\n",
"                end\n",
"                if(k<l)\n",
"                    if(k==j & l==i)\n",
"                        C(k,l)=s(i,j)\n",
"                    else\n",
"                        C(k,l)=0\n",
"                    end\n",
"                end\n",
"            end\n",
"        end\n",
"        printf('\n\n Iteration %i',m)\n",
"        m=m+1\n",
"        disp(C,'C=');\n",
"        R=C*R;\n",
"        Q=Q*C';\n",
"        disp(Q,'Q=',R,'R=')\n",
"    end\n",
"end\n",
"disp(Q*R,'Q*R=A=')            //verification"
   ]
   }
,
{
		   "cell_type": "markdown",
		   "metadata": {},
		   "source": [
			"## Example 9.8: Recursive_Least_Square_Method.sce"
		   ]
		  },
  {
"cell_type": "code",
	   "execution_count": null,
	   "metadata": {
	    "collapsed": true
	   },
	   "outputs": [],
"source": [
"//Example 9.8\n",
"//Recursive Least-Square Method\n",
"//Page no. 308\n",
"clc;clear;close;\n",
"\n",
"A0=[3,0;0,3;3,3];\n",
"B0=[2;2;2];\n",
"A1=[6,3];B1=[6];\n",
"A0T=A0';\n",
"G0=A0T*A0;\n",
"disp(G0,'G0=')\n",
"G0_1=inv(G0);\n",
"disp(G0_1,'Inverse of G0=')\n",
"X0=G0_1*A0T*B0;\n",
"disp(X0,'X0=')\n",
"\n",
"//by recursive least square algorithm\n",
"G1=G0+A1'*A1;\n",
"disp(G1,'G1=');\n",
"G1_1=inv(G1);\n",
"disp(G1_1,'Inverse of G1')\n",
"X1=X0+G1_1*A1'*(B1-A1*X0);\n",
"disp(X1,'X1=')\n",
"\n",
"//verification\n",
"A=[3,0;0,3;3,3;6,3];\n",
"B=[2;2;2;6];\n",
"AT=A';\n",
"G=AT*A;\n",
"disp(G,'G=')\n",
"G_1=inv(G);\n",
"disp(G_1,'Inverse of G=')\n",
"X=G_1*AT*B;\n",
"disp(X,'X=')\n",
"disp('Thus X and X1 are Same')"
   ]
   }
],
"metadata": {
		  "kernelspec": {
		   "display_name": "Scilab",
		   "language": "scilab",
		   "name": "scilab"
		  },
		  "language_info": {
		   "file_extension": ".sce",
		   "help_links": [
			{
			 "text": "MetaKernel Magics",
			 "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
			}
		   ],
		   "mimetype": "text/x-octave",
		   "name": "scilab",
		   "version": "0.7.1"
		  }
		 },
		 "nbformat": 4,
		 "nbformat_minor": 0
}