summaryrefslogtreecommitdiff
path: root/Numerical_Methods_Principles_by_Analysis/9-Linear_Least_Squares_Problem.ipynb
diff options
context:
space:
mode:
Diffstat (limited to 'Numerical_Methods_Principles_by_Analysis/9-Linear_Least_Squares_Problem.ipynb')
-rw-r--r--Numerical_Methods_Principles_by_Analysis/9-Linear_Least_Squares_Problem.ipynb479
1 files changed, 479 insertions, 0 deletions
diff --git a/Numerical_Methods_Principles_by_Analysis/9-Linear_Least_Squares_Problem.ipynb b/Numerical_Methods_Principles_by_Analysis/9-Linear_Least_Squares_Problem.ipynb
new file mode 100644
index 0000000..0e7132d
--- /dev/null
+++ b/Numerical_Methods_Principles_by_Analysis/9-Linear_Least_Squares_Problem.ipynb
@@ -0,0 +1,479 @@
+{
+"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
+}
generated by cgit (git 2.34.1) at 2025-07-13 01:30:59 +0000