Introduction to Matlab Programming with Applications

Size: px

Start display at page:

Download "Introduction to Matlab Programming with Applications"

综晖帛
4 years ago
Views:

1 1 >> Lecture 7 2 >> 3 >> -- Matrix Computation 4 >> Zheng-Liang Lu 333 / 410

2 Vectors ˆ Let R be the set of all real numbers. ˆ R m 1 denotes the vector space of all m-by-1 column vectors: u = (u i ) = u 1. u m ˆ Similarly, the row vector v R 1 n is. (1) v = (v i ) = [ v 1 v n ]. (2) Zheng-Liang Lu 334 / 410

3 Matrices ˆ R m n denotes the vector space of all m-by-n real matrices A: a 11 a 1n A = (a ij ) = a m1 a mn ˆ Recall that we use the subscripts to reference a particular element in a matrix. ˆ The definition for complex matrices is similar. Zheng-Liang Lu 335 / 410

4 Transposition ˆ transpose(a) returns the transpose of A. 1 >> A = magic(4) % a special matrix provided by Matlab 2 >> transpose(a) ˆ You can also transpose A by 1 >> A' ˆ Note that if A C m n, then A returns the transposition and complex conjugate 1 of A. 1 z = a + bi has a complex conjugate z = a bi. Zheng-Liang Lu 336 / 410

5 Arithmetic Operations ˆ Let a ij and b ij be the elements of the matrices A and B R m n for 1 i m and 1 j n, respectively. ˆ Then C = A ± B can be calculated by c ij = a ij ± b ij. (Try.) Zheng-Liang Lu 337 / 410

6 ˆ Let u, v R m 1. Inner Product 2 ˆ Then the inner product, denoted by u v, is calculated by u T v = [u 1 u m ] where T is the transposition operator. v 1. v m, ˆ Inner product is also called projection project for emphasizing the geometric significance. 2 Akaa dot product and scalar product. Zheng-Liang Lu 338 / 410

7 Example 1 clear; clc; 2 3 x = [1; 2; 3]; 4 y = [4; 5; 6]; 5 z = 0; 6 for i = 1 : 3 7 z = z + x(i) * y(i); 8 end 9 z % by definition; using a for loop 10 x' * y % normal way 11 dot(x, y) % using built-in function Zheng-Liang Lu 339 / 410

8 Exercise ˆ Let u and v be any two vectors in R m 1. ˆ Write a program that returns the angle between u and v. Input: u, v Output: θ in radian ˆ You may use norm(u) to calculate the length of u. 3 ˆ Alternatively, use subspace(u, v). 3 See Zheng-Liang Lu 340 / 410

9 Generalization of Inner Product ˆ For simplicity, consider x R. ˆ Let f (x) be a real-valued function. ˆ Let g(x) be a basis function. 4 ˆ Then we can define the inner product of f and g on [a, b] by f, g = b a f (x)g(x)dx. ˆ There exist plenty of basis functions which can be used to represent the functions. 5 4 See 5 See and Zheng-Liang Lu 341 / 410

10 ˆ For example, Fourier transform is widely used in engineering and science. ˆ Fourier integral 6 is defined as F (ω) = f (t)e iωt dt where f (t) is a square-integrable function. ˆ The Fast Fourier transform (FFT) algorithm computes the discrete Fourier transform (DFT) in O(n log n) time. 7,8 6 See 7 Cooley and Tukey (1965). 8 See Zheng-Liang Lu 342 / 410

11 Cross Product 10 ˆ cross(u, v) returns the cross product of the vectors x and y of length 3. 9 ˆ For example, 1 >> x = [1; 0; 0]; 2 >> y = [0; 1; 0]; 3 >> z = cross(x, y) % built-in function 4 5 z = Actually, only in 3- and 7-dimensional Euclidean spaces. 10 For example, angular momentum, Lorentz force, and Poynting vector. Zheng-Liang Lu 343 / 410

12 Zheng-Liang Lu 344 / 410

13 Matrix Multiplication ˆ Let A R m q and B R q n. ˆ Then C = A B is given by c ij = ˆ For example, q a ik b kj. (3) k=1 Zheng-Liang Lu 345 / 410

14 Example 1 clear; clc; 2 3 A = randi(10, 5, 4); 4 B = randi(10, 4, 3); 5 C = zeros(size(a, 1), size(b, 2)); 6 for i = 1 : size(a, 1) 7 for j = 1 : size(b, 2) 8 for k = 1 : size(a, 2) 9 C(i, j) = C(i, j) + A(i, k) * B(k, j); 10 end 11 end 12 end 13 C ˆ Time complexity: O(n 3 ) ˆ Strassen (1969): O(n ) ˆ Coppersmith and Winograd (2010): O(n ) Zheng-Liang Lu 346 / 410

15 Matrix Exponentiation ˆ Raising a matrix to a power is equivalent to repeatedly multiplying the matrix by itself. ˆ For example, A 2 = A A. ˆ It implies that A should be square. (Why?) ˆ The matrix exponential 11 is a matrix function on square matrices analogous to the ordinary exponential function. ˆ More explicitly, e A A n = n!. ˆ However, raising a matrix to a matrix power, that is, A B, is undefined. n=0 11 See matrix exponentials and Pauli matrices. Zheng-Liang Lu 347 / 410

16 Example ˆ Be aware that exp(a) and exp(1) A are different for any square matrix A. 1 clear; clc; 2 3 A = [0 1; 1 0]; 4 a = exp(1) ˆ A 5 b = eye(size(a)); 6 for n = 1 : 10 7 b = b + X ˆ n / factorial(n); 8 end 9 b 10 c = exp(a) % different? Zheng-Liang Lu 348 / 410

17 1 a = b = c = ˆ Use expm(a) and logm(a) for square A. Zheng-Liang Lu 349 / 410

18 Determinants ˆ det(a) denotes the determinant of the square matrix A. [ ] a b ˆ Recall that if A =, then det(a) = ad bc. c d ˆ The function det(a) returns the determinant of A. ˆ For example, 1 >> A = magic(3); 2 >> det(a) ˆ However, the general formula 12 is somehow complicated. 12 See Zheng-Liang Lu 350 / 410

19 Recursive Algorithm for det 1 function y = mydet(a) 2 [r, c] = size(a); 3 if r == c % check r == c? 4 if c == 1 5 y = A; 6 elseif c == 2 7 y = A(1, 1) * A(2, 2) - A(1, 2) *... A(2, 1); 8 else 9 y = A(1, 1) * mydet(a(2 : c, 2 : c)); 10 for i = 2 : c y = y + A(1, i) * (-1) ˆ (i + 1)... * mydet(a(2 : c, [1 : i - 1, i : c])); 12 end 13 y = y + A(1, c) * (-1) ˆ (c + 1) *... mydet(a(2 : c, 1 : c - 1)); Zheng-Liang Lu 351 / 410

20 14 end 15 else 16 error('should be a square matrix.') 17 end 18 end ˆ Clearly, there are n! terms in the calculation. ˆ This algorithm runs in O(2 n n 2 ) time! ˆ An efficient calculation of determinants can be done by LU decomposition which runs in O(n 3 ) time See LU decomposition, QR factorization, and Cholesky factorization. Zheng-Liang Lu 352 / 410

21 Inverse Matrices ˆ A R n n is invertible if there exists B R n n such that A B = B A = I n, where I n denotes the n-by-n identity matrix. ˆ You can use eye(n) to generate an identity matrix I n. ˆ A is invertible if and only if det(a) 0. ˆ The function inv(a) returns the inverse of A. Zheng-Liang Lu 353 / 410

22 ˆ For example, 1 >> A = pascal(4); % Pascal matrix 2 >> B = inv(a) ˆ However, inv(a) may return a result even if A is ill-conditioned. 14 ˆ For example, det = Try rcond. Zheng-Liang Lu 354 / 410

23 First-Order Approximation 16 ˆ Let f (x) be a nonlinear function which is infinitely differentiable at x 0. ˆ By Taylor s expansion 15, we have f (x) f (x 0 ) + f (x 0 )(x x 0 ) + O( x 2 ), where x 2 is the higher-order term, which can be neglected as x x 0. ˆ So a nonlinear function can be reduced by its first-order approximation. 15 See 16 For example, the Newtonian equation is only a low-speed approximation to Einstein s total energy, given by E = m 0c 1 (v/c) 2, where m0 is the rest mass and 2 v is the velocity relative to the inertial coordinate. By applying the first-order approximation, E m 0c mv 2. 2 Zheng-Liang Lu 355 / 410

24 System of Linear Equations ˆ A system of linear equations is a set of linear equations involving the same set of variables. ˆ For example, nodal analysis by Kirchhoff s Laws. Zheng-Liang Lu 356 / 410

25 ˆ A general system of m linear equations with n unknowns can be written as a 11 x 1 +a 12 x 2 +a 1n x n = b 1 a 21 x 1 +a 22 x 2 +a 2n x n = b =. a m1 x 1 +a m2 x 2 +a mn x n = b m where x 1,..., x n are unknowns, a 11,..., a mn are the coefficients of the system, and b 1,..., b m are the constant terms. Zheng-Liang Lu 357 / 410

26 ˆ Hence we can rewrite the system of linear equations as a matrix equation, given by Ax = b. where a 11 a 12 a 1n a 21 a 22 a 2n A =......, a m1 a m2 a mn x = x 1. x n, and b = b 1. b m. Zheng-Liang Lu 358 / 410

27 Solving General System of Linear Equations ˆ Let x be the column vector with n independent variables and m constraints. 17 ˆ If m = n, then there exists the unique solution. 18 ˆ If m > n, then there is no exact solution. ˆ Fortunately, we can find a least-squares error solution such that Ax b 2 is minimal. (See the next page.) ˆ If m < n, then there are infinitely many solutions. ˆ We can calculate the inverse by simply using the left matrix divide operator (\) or mldivide like this: x = A\b. 17 Assume that they are linearly independent. 18 Equivalently, rank(a) = rank([a, b]). Also see Cramer s rule. Zheng-Liang Lu 359 / 410

28 Unique Solution (m = n) ˆ For example, 3x +2y z = 1 x y +2z = 1 2x +y 2z = 0 1 >> A = [3 2-1; 1-1 2; ]; 2 >> b = [1; -1; 0]; 3 >> x = A \ b Zheng-Liang Lu 360 / 410

29 Overdetermined System (m > n) ˆ For example, 2x y = 2 x 2y = 2 x +y = 1 1 >> A=[2-1; 1-2; 1 1]; 2 >> b=[2; -2; 1]; 3 >> x = A \ b Zheng-Liang Lu 361 / 410

30 ˆ For example, Underdetermined System (m < n) { x +2y +3z = 7 4x +5y +6z = 8 1 >> A = [1 2 3; 4 5 6]; 2 >> b = [7; 8]; 3 >> x = A \ b % (Why?) ˆ Note that this solution is a basic solution, one of infinitely many. ˆ How to find the directional vector? Zheng-Liang Lu 362 / 410

31 Gaussian Elimination ˆ Recall the procedure of Gaussian Elimination in high school. ˆ Now we proceed to write a program which solves the following simultaneous equations: 3x +2y z = 1 x y +2z = 1 2x +y 2z = 0 ˆ Then we have x = 1, y = 2, and z = 2. Zheng-Liang Lu 363 / 410

32 ˆ Suppose det(a) 0. ˆ Form an upper triangular matrix 1 ā 12 ā 1n 0 1 ā 2n Ā =.. 1. with b = where ā ij s and b i s are the values after math. b 1 b 2. b n, ˆ Use a backward substitution to determine the solution vector x by where i = 1, 2,, n. x i = b i n j=i+1 ā ij x j, Zheng-Liang Lu 364 / 410

33 Solution 1 clear; clc; 2 3 A = [3 2-1; 1-1 2; ]; 4 b = [1; -1; 0]; 5 A \ b % check the answer 6 7 if det(a) ~= 0 8 for i = 1 : 3 9 for j = i : 3 10 % cannot be interchanged % 11 b(j) = b(j) / A(j, i); 12 A(j, :) = A(j, :) / A(j, i); 13 % % % % % % % % % % % % % % 14 end 15 for j = i + 1 : 3 16 A(j, :) = A(j, :) - A(i, :); 17 b(j) = b(j) - b(i); Zheng-Liang Lu 365 / 410

34 18 end 19 end 20 x = zeros(3, 1); 21 for i = 3 : -1 : 1 22 x(i) = b(i); 23 for j = i + 1 : 1 : 3 24 x(i) = x(i) - A(i, j) * x(j); 25 end 26 end 27 else 28 disp('no unique solution.'); 29 end 30 x Zheng-Liang Lu 366 / 410

35 Exercise ˆ Write a program which solves a general system of linear equations. ˆ The function rank(a) provides an estimate of the number of linearly independent rows or columns of A. 19 ˆ Check if rank(a) = rank([a, b]). ˆ If so, then there is at least one solution. ˆ If not, then there is no solution. ˆ The function rref([a, b]) produces the reduced row echelon form of A. 19 rank(a) min{r, c} where r and c are the numbers of rows and columns. Zheng-Liang Lu 367 / 410

36 Zheng-Liang Lu 368 / 410

37 Solution 1 function y = linearsolver(a, b) 2 3 if rank(a) == rank([a, b]) % argumented matrix 4 if rank(a) == size(a, 2); 5 disp('exact one solution.') 6 x = A \ b 7 else 8 disp('infinite numbers of solutions.') 9 rref([a b]) 10 end 11 else 12 disp('there is no solution. (Only least... square solutions.)') 13 end Zheng-Liang Lu 369 / 410

38 Example: 2D Laplace s Equation for Electrostatics ˆ Laplace s equation 20 is one of 2nd-order partial differential equations (PDEs). 21 ˆ Let Φ(x, y) be an electrical potential, which is a function of x, y R. ˆ Consider 2 Φ(x, y) = 0, where 2 = is the Laplace operator. x 2 y 2 ˆ Solving Laplace s equation in practical applications often requires numerical methods. 20 Pierre-Simon Laplace ( ). 21 See Zheng-Liang Lu 370 / 410

39 Rectangular Trough Zheng-Liang Lu 371 / 410

40 Extremely Simple Assumption ˆ First, we can partition the region into many subregions by a proper mesh generation. ˆ If Φ(x, y) satisfies the Laplace s equation, then Φ(x, y) can be approximated by Φ(x, y) Φ(x + ɛ, y) + Φ(x ɛ, y) + Φ(x, y + ɛ) + Φ(x, y ɛ), 4 where ɛ is a small distance compared with the system size. Zheng-Liang Lu 372 / 410

41 V5 V10 V15 V20 V V4 V9 V14 V19 V V3 V8 V13 V18 V V2 V7 V12 V17 V V1 V V11 V16 V21 Zheng-Liang Lu 373 / 410

42 Reformulation ˆ Consider the boundary condition: ˆ V 1 = V 2 = = V 4 = 0 ˆ V 21 = V 22 = = V 24 = 0 ˆ V 1 = V 6 = = V 16 = 0 ˆ V 5 = V 10 = = V 25 = 100 ˆ Now define x = [ V 7 V 8 V 9 V 12 V 13 V 14 V 17 V 18 V 19 ] T where T is the transposition operator. Zheng-Liang Lu 374 / 410

43 ˆ Then we form Ax = b where A = and b = [ ] T. ˆ As you can see that V 7 = V 17, V 8 = V 18 and V 9 = V 19 due to the spatial symmetry, the dimension of A can be reduced to 6! (Try.) Zheng-Liang Lu 375 / 410

44 1 clear; clc; close all; 2 3 a = 1; b = 1; n = 5; V0 = 100; 4 5 x = linspace(0, a, 5); 6 y = linspace(0, b, 5); 7 [X Y] = meshgrid(x, y); 8 9 figure; hold on; grid on; 10 plot(x, Y, 'k.', 'markersize', 24); 11 for i = 1 : length(x) 12 for j = 1 : length(y) 13 text(x(n * (i - 1) + j), Y(n * (i - 1) +... j) , sprintf('v%d', n * (i - 1)... + j)); 14 end 15 end % boundary condition Zheng-Liang Lu 376 / 410

45 18 phi = zeros(1, length(x) * length(y)); 19 phi(5 : 5 : 25) = 100; A = [ ; ; ; ; ; ]; 27 bb = [0; 0; 100; 0; 0; 100]; % inverse of the matrix 30 v = A \ bb; % generate the solution matrix 33 phi([7 8 9]) = v(1 : 3); 34 phi([ ]) = phi([7 8 9]); 35 phi([ ]) = v(4 : 6); phi = reshape(phi, 5, 5); 38 for i = 1 : length(y) Zheng-Liang Lu 377 / 410

46 39 for j = 1 : length(x) 40 h = text(x(n * (i - 1) + j), Y(n * (i ) + j) , sprintf('%7.4f',... phi(j, i))); 41 set(h, 'color', 'r'); 42 end 43 end figure; hold on; grid on; 46 contour(x, Y, phi); colorbar; ˆ This is a toy example for numerical methods. ˆ You may consider Finite Difference Method (FDM) and Finite Element Method (FEM), both widely used in commercial simulation softwares! 22 ˆ Besides, the mesh generation is also important for numerical methods Read 23 See Zheng-Liang Lu 378 / 410

47 Zheng-Liang Lu 379 / 410

48 Method of Least Squares ˆ The first clear and concise exposition of the method of least squares was published by Legendre in ˆ In 1809, Gauss published his method of calculating the orbits of celestial bodies. ˆ The method of least squares is a standard approach to the approximate solution of overdetermined systems, that is, sets of equations in which there are more equations than unknowns. 24 ˆ To obtain the coefficient estimates, the least-squares method minimizes the summed square of residuals. 24 Aka degrees of freedom. Zheng-Liang Lu 380 / 410

49 ˆ Let {ŷ i } n i=1 be the observed response values and {y i} n i=1 be the fitted response values. ˆ Let ɛ i = ŷ i y i be the residual for i = 1,..., n. ˆ Then the sum of square error estimates associated with the data is given by n S = ɛ 2 i. i=1 Zheng-Liang Lu 381 / 410

50 Illustration Zheng-Liang Lu 382 / 410

51 Linear Least Squares ˆ In the sense of linear least squares, a linear model is said to be an equation which is linear in the coefficients. ˆ Now we choose a linear equation, y = ax + b, where a and b are to be determined. ˆ So ɛ i = (ax i + b) ŷ i and then S = n ((ax i + b) ŷ i ) 2. i=1 ˆ The coefficient a and b can be determined by differentiating S with respect to each parameter, and setting the result equal to zero. (Why?) Zheng-Liang Lu 383 / 410

52 ˆ More explicitly, S n a = 2 x i (y i (ax i + b)) = 0, i=1 S n b = 2 (y i (ax i + b)) = 0. i=1 ˆ So the aforesaid equations are reorganized as a n xi 2 i=1 a n + b x i = i=1 n x i + nb = n x i y i, i=1 i=1 i=1 n y i. Zheng-Liang Lu 384 / 410

53 ˆ In form of matrices, [ n i=1 x i 2 n i=1 x i n i=1 x i n ˆ So we have ] [ a b ] = [ n i=1 x iy i n i=1 y i ]. a = n n i=1 x iy i n i=1 x n i i=1 y i n n i=1 x i 2 ( n i=1 x i) 2 = cov(x, y) cov(x), where cov(x, y) denotes the covariance between x = {x i } n i=1 and y = {y i } n i=1. ˆ Then we have b = 1 n n ( y i a i=1 n x i ). i=1 Zheng-Liang Lu 385 / 410

54 Example: Circle Fitting ˆ Consider a set of data points surrounding some center. ˆ Now the coordinates of the circle center and also the radius are desired. ˆ This needs to estimate 3 unknowns: (x c, y c ) and r > 0. ˆ Recall that a circle equation is (x x c ) 2 + (y y c ) 2 = r 2. ˆ The above equation can be equivalent to 2xx c + 2yy c + z = x 2 + y 2, where z = r 2 x 2 c + y 2 c. Zheng-Liang Lu 386 / 410

55 ˆ For a set of data points (x i, y i ), i = 1, 2, 3,..., N, this rearranged equation can be written in matrix form Aw = b, where 2x 1 2y 1 1 x c x1 2 A =....., w = y c + y 2 1, b =.. 2x N 2y N 1 z xn 2 + y N 2 Zheng-Liang Lu 387 / 410

56 1 clear; clc; close all; 2 3 N = 100; 4 theta = 2 * pi * rand(1, N); 5 xcc = 5; 6 ycc = 3; 7 rcc = 10; 8 9 x = xcc + rcc * cos(theta) + randn(1, N) * 0.5; 10 y = ycc + rcc * sin(theta) + randn(1, N) * 0.5; xt = x - mean(x); 13 yt = y - mean(y); 14 distance = sqrt(xt.ˆ 2 + yt.ˆ 2) 15 maxr = max(distance); xt = xt / maxr; 18 yt = yt / maxr; 19 distance = distance / maxr; Zheng-Liang Lu 388 / 410

57 20 21 A = [2 * xt', 2 * yt', ones(n, 1)]; 22 b = (distance.ˆ 2)'; % v = [xc; yc; z] 25 v = A \ b 26 r = sqrt(v(3) + v(1) ˆ 2 + v(2) ˆ 2) * maxr 27 xc = v(1) * maxr + mean(x) 28 yc = v(2) * maxr + mean(y) figure; plot(x, y, 'o'); 31 hold on; grid on; axis equal; theta = linspace(0, 2 * pi, 100); 34 x = xc + r * cos(theta ); 35 y = yc + r * sin(theta ); 36 plot(x, y, 'r-'); Zheng-Liang Lu 389 / 410

58 Zheng-Liang Lu 390 / 410

59 Polynomials 25 ˆ In fact, all polynomials of n-th order with addition and multiplication to scalars form a vector space, denoted by P n. ˆ In general, f (x) is said to be a polynomial of n-order provided that f (x) = a n x n + a n 1 x n a 0, where a n 0. ˆ It is convenient to express a polynomial by a coefficient vector (a n, a n 1,..., a 0 ), where the elements are the coefficients of the polynomial in descending order. 25 Weierstrass approximation theorem states that every continuous function defined on a closed interval [a, b] can be uniformly approximated as closely as desired by a polynomial function. See Zheng-Liang Lu 391 / 410

60 Arithmetic Operations ˆ P 1 + P 2 returns the addition of two polynomials. ˆ P 1 P 2 returns the subtraction of two polynomials. ˆ The function conv(p 1, P 2 ) returns the resulting coefficient vector for multiplication of the two polynomials P 1 and P ˆ The function [Q, R] = deconv(b, A) deconvolves vector A out of vector B. ˆ Equivalently, B = conv(a, Q) + R. ˆ This is so-called Euclidean division algorithm. ˆ The function polyval(p, X ) returns the values of a polynomial P evaluated at x X. 26 See Convolution. Zheng-Liang Lu 392 / 410

61 1 clear; clc; 2 3 p1 = [ ]; 4 p2 = [ ]; 5 %%% addition 6 p3 = p1 + p2 7 %%% substraction 8 p4 = p1 - p2 9 %%% multiplcaition 10 p5 = conv(p1, p2) 11 %%% division: q is quotient and r is remainder 12 [q, r] = deconv(p1, p2) 13 x = -1 : 0.1 : 1; 14 plot(x, polyval(p1, x), 'o', x, polyval(p2, x),... '*', x, polyval(p5, x), 'd'); 15 grid on; legend('p1', 'p2', 'conv(p1, p2)'); Zheng-Liang Lu 393 / 410

62 30 p1 p2 conv(p1,p2) Zheng-Liang Lu 394 / 410

63 Roots Finding ˆ The function roots(p) returns a vector whose elements are all roots of the polynomial P. 27 ˆ For example, 1 clear; clc; 2 3 p = [1, 3, 1, 5, -1]; 4 r = roots(p) 5 x = -4 : 0.1 : 1; 6 plot(x, polyval(p, x), '--'); hold on; grid on; 7 for i = 1 : length(r) 8 if isreal(r(i)) == 1 9 plot(r, polyval(p, r(i)), 'ro'); 10 end 11 end 12 polyval(p, r) 27 See Zheng-Liang Lu 395 / 410

64 1 >> r = i i >> ans = e-013 * i i 15 0 ˆ Why not exactly zero? Zheng-Liang Lu 396 / 410

65 Zheng-Liang Lu 397 / 410

66 Exercise: Internal Rate of Return (IRR) ˆ Given a collection of pairs (time, cash flow) involved in a project, the IRR is a rate of return when the net present value is zero. ˆ Explicitly, the IRR can be calculated by solving N n=0 C n (1 + r) n = 0, where C n is the cash flow at time n. ˆ For example, consider an investment may be given by the sequence of cash flows: C 0 = , C 1 = 36200, C 2 = 54800, C 3 = ˆ Then the IRR is 5.96%. Zheng-Liang Lu 398 / 410

67 Forming Polynomials ˆ The function poly(v ), where V is a vector, returns a vector whose elements are the coefficients of the polynomial whose roots are the elements of V. ˆ Simply put, the function roots and poly are inverse functions of each other. Zheng-Liang Lu 399 / 410

68 Example 1 clear; clc; 2 3 v = [0.5 sqrt(2) 3]; 4 y = 1; 5 for i = 1 : 3 6 y = conv(y, [1 -v(i)]); 7 end 8 y 9 10 poly(v) Zheng-Liang Lu 400 / 410

69 Integral and Derivative of Polynomials ˆ The function polyder(p) returns the derivative of the polynomial whose coefficients are the elements of vector P in descending powers. ˆ The function polyint(p, K) returns a polynomial representing the integral of polynomial P, using a scalar constant of integration K. 1 clear; clc; 2 3 p = [ ]; 4 p der = polyder(p) 5 p int = polyint(p, 0) % assume K = 0 Zheng-Liang Lu 401 / 410

70 Exercise ˆ Consider f (x) = 4x 3 + 3x 2 + 2x + 1 for x R. ˆ Determine the coefficients of its derivative f and integration F (x) = x 0 f (t)dt. ˆ Do not use the built-in functions. Zheng-Liang Lu 402 / 410

71 1 clear; clc; 2 3 p = [ ]; 4 K = 0; 5 q1 = zeros(1, length(p)); 6 for i = 2 : length(p) q1(i) = p(i - 1) * (length(p) - (i - 1)); 8 end 9 q q2 = zeros(1, length(p) + 1); 12 q2(length(q2)) = K; 13 for i = 1 : length(p) 14 q2(i) = 1 / (length(p) - i + 1) * p(i); 15 end 16 q2 Zheng-Liang Lu 403 / 410

72 Curve Fitting by Polynomials ˆ The function polyfit(x, y, n) returns the coefficients for a polynomial p(x) of degree n that is a best fit (in a least-squares sense) for the data in y. Zheng-Liang Lu 404 / 410

73 Example 1 clear; clc; close all; 2 3 x = linspace(0, 1, 10); 4 y = cos(rand(1, length(x)) * pi / 2) + x.ˆ 2; 5 figure; hold on; grid on; plot(x, y, 'o'); 6 7 color = 'rgbck'; 8 x = linspace(0, 1, 100); 9 for i = 1 : 5 10 p = polyfit(x, y, i); 11 plot(x, polyval(p, x ), color(i)); 12 end 13 p A = [x'.ˆ 5, x'.ˆ 4, x'.ˆ 3, x'.ˆ 2, x'.ˆ 1,... ones(10, 1)]; 16 b = y'; 17 pp = A \ b Zheng-Liang Lu 405 / 410

74 Overfitting Zheng-Liang Lu 406 / 410

75 Occam s Razor Entities must not be multiplied beyond necessity. Duns Scotus ˆ In science, Occam s razor is used as a heuristic to guide scientists in developing theoretical models rather than as an arbiter between published models. ˆ Among competing hypotheses, the one with the fewest assumptions should be selected. ˆ For example, Runge s phenomenon is a problem of oscillation at the edges of an interval that occurs when using polynomial interpolation with polynomials of high degree over a set of equispaced interpolation points See s_phenomenon. Zheng-Liang Lu 407 / 410

76 Eigenvalues and Eigenvectors 29 ˆ Let A be a square matrix. ˆ Then v is an eigenvector associated with the eigenvalue λ if ˆ Equivalently, Av = λv. (A λi )v = 0. ˆ For nontrivial vectors v, det(a λi ) = 0. ˆ The above equation is the so-called characteristic polynomial, whose roots are actually eigenvalues! ˆ Use eig(a) to derive the eigenvalues associated with eigenvectors for the matrix A. 29 See eigenvectors#applications. Zheng-Liang Lu 408 / 410

77 Singular Value Decomposition (SVD) 30 ˆ Let A m n be a matrix. ˆ Then σ is called one singular value associated with the singular vectors u R m 1 and v R n 1 for A provided that { Av = σu, A T u = σv. ˆ We further have { AV = UΣ, A T U = V Σ, where U and V are both unitary, and the diagonal terms in Σ are σ s, 0 s in off-diagonal terms. ˆ You may use the built-in function svd. 30 See Zheng-Liang Lu 409 / 410

78 Example: Low-rank Approximation for Image Compression ˆ This idea originates from Principal Component Analysis (PCA). 31 ˆ Use svd to calculate the principal components of the input image. ˆ Then we can have an image extremely similar to the origin one, but with a smaller image size by keeping the vectors associated with a few first largest of principal components. 31 See PCA-Tutorial-Intuition_jp.pdf and Zheng-Liang Lu 410 / 410

Introduction to Matlab Programming with Applications

Introduction to Matlab Programming with Applications 1 >> Lecture 7 2 >> 3 >> -- Matrix Computation 4 >> Zheng-Liang Lu 297 / 343 Vectors ˆ Let R be the set of all real numbers. ˆ R m 1 denotes the vector space of all m-by-1 column vectors: u = (u i ) =