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 297 / 343

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 298 / 343

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 299 / 343

4 Transposition 1 >> A = [1 i]; 2 >> A' % Hermitian operator; see any textbook for... linear algebra 3 4 ans = i i >> A.' % transposition of A 11 ans = i i Zheng-Liang Lu 300 / 343

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. ˆ Then C = A ± B can be calculated by c ij = a ij ± b ij. (Try.) Zheng-Liang Lu 301 / 343

6 ˆ Let u, v R m 1. Inner Product 1 ˆ Then the inner product, denoted by u v, is calculated by u v = u v = [u 1 u m ] v 1. v m. ˆ Inner product is also called projection for emphasizing the geometric significance. 1 Akaa dot product and scalar product. Zheng-Liang Lu 302 / 343

7 Example 1 clear; clc; 2 3 u = [1; 2; 3]; 4 v = [4; 5; 6]; 5 u' * v % normal way 6 dot(u, v) % using the built-in function Zheng-Liang Lu 303 / 343

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. 2 ˆ Alternatively, use subspace(u, v). 2 See Zheng-Liang Lu 304 / 343

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. 3 ˆ 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. 4 3 See 4 See and Zheng-Liang Lu 305 / 343

10 ˆ For example, Fourier transform is widely used in engineering and science. ˆ Fourier integral 5 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. 6,7 5 See 6 Cooley and Tukey (1965). 7 See Zheng-Liang Lu 306 / 343

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

12 Zheng-Liang Lu 308 / 343

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 309 / 343

14 Example 1 clear; clc; 2 3 A = randi(10, 5, 4); % 5-by-4 4 B = randi(10, 4, 3); % 4-by-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 % display C ˆ Time complexity: O(n 3 ) ˆ Strassen (1969): O(n ) Zheng-Liang Lu 310 / 343

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 10 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 10 See matrix exponentials and Pauli matrices. Zheng-Liang Lu 311 / 343

16 Example 1 clear; clc; 2 3 A = [0 1; 1 0]; 4 5 trial1 = exp(a) 6 trial2 = eye(size(a)); 7 for n = 1 : 10 8 trial2 = trial2 + X ˆ n / factorial(n); 9 end 10 trial3 = exp(1) ˆ A % equivalent to expm(a) Zheng-Liang Lu 312 / 343

17 1 trial1 = trial2 = trial3 = Zheng-Liang Lu 313 / 343

18 Determinants ˆ det(a) denotes the determinant of A. [ ] a b ˆ Given A =, then c d det(a) = ad bc. ˆ The function det(a) returns the determinant of A. ˆ However, the general formula for det(a) 11 is complicated. 11 See Zheng-Liang Lu 314 / 343

19 Naive Algorithm 1 function y = mydet(a) 2 [r, ~] = size(a); 3 4 if r == 1 5 y = A; 6 elseif r == 2 7 y = A(1, 1) * A(2, 2) - A(1, 2) * A(2, 1); 8 else 9 y = 0; 10 for i = 1 : r 11 adjoint = A(2 : r, [1 : i - 1, i : r]); 12 y = y + A(1, i) * (-1) ˆ (i + 1) *... mydet(adjoint); 13 end 14 end 15 end Zheng-Liang Lu 315 / 343

20 ˆ It needs n! terms in the sum of products, so this algorithm runs in O(n!) time! ˆ You could compare this algorithm to det for a 10-by-10 matrix. ˆ In practice, an efficient calculation of determinants can be done in O(n 3 ) time, say by LU decomposition See Zheng-Liang Lu 316 / 343

21 Inverse Matrices 13 ˆ 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. 13 See invertible_matrix_theorem. Zheng-Liang Lu 317 / 343

22 ˆ For example, 1 >> A = pascal(4); % Pascal matrix 2 >> inv(a) % Try it! ˆ However, inv(a) may return a result even if A is ill-conditioned. 14 ˆ For example, det = You may refer to the condition number of a function with respect to an argument. This number measures how much the output value of the function can change for a small change in the input argument. Try rcond. Zheng-Liang Lu 318 / 343

23 First-Order Approximation ˆ 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( 2 x), where 2 x = (x x 0 ) 2 is the higher-order term, which can be neglected as x 0. ˆ Then we have f (x) f (x 0 )x + k, with k = f (x 0 ) f (x 0 )x 0 is a constant. 15 See Zheng-Liang Lu 319 / 343

24 A Profound Concept: Local Linearization ˆ For example, we barely feel like the curvature of the ground; however, we watch Earth on the moon and agree that Earth is a sphere. ˆ For another example, the Newtonian equation is only a low-speed approximation to Einstein s total energy. ˆ Let m 0 be the rest mass and v be the velocity relative to the inertial coordinate. ˆ Then the total energy is E = m 0 c 2 1 (v/c) 2. ˆ By applying the first-order approximation, E m 0 c mv 2. Zheng-Liang Lu 320 / 343

25 System of Linear Equations ˆ A system of linear equations is a collection of two or more linear equations involving the same set of variables. ˆ For example, Kirchhoff s Laws. Zheng-Liang Lu 321 / 343

26 General Form ˆ 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 322 / 343

27 Matrix Equation ˆ Hence we can rewrite the aforesaid equations as follows: 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 323 / 343

28 Solutions to System of Linear Equations 17 ˆ Let n be the number of unknowns and m be the number of constraints. 16 ˆ If m = n, then there exists a single unique solution. ˆ If m > n, then there is no solution. ˆ Fortunately, we can find a least-squares error solution such that Ax b 2 is minimal. ˆ If m < n, then there are infinitely many solutions. ˆ We can calculate solutions of these three kinds by x = A \ b. 16 Here we assume that these m row vectors from the constraints are linearly independent. See 17 See Zheng-Liang Lu 324 / 343

29 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 325 / 343

30 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 326 / 343

31 ˆ 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 ˆ Note that this solution is a basic solution, one of infinitely many. ˆ How to find the directional vector? Zheng-Liang Lu 327 / 343

32 Gaussian Elimination Algorithm 18 ˆ First we consider the linear system is represented as an augmented matrix [ A b ]. ˆ We then transform A into an upper triangular matrix 1 ā 12 ā 1n b1 0 1 ā 2n b 2 [ Ā b ] = bn where ā ij s and b i s are the resulting values after elementary row operations. ˆ This matrix is said to be in reduced row echelon form. 18 See Zheng-Liang Lu 328 / 343

33 ˆ By this form, we can have the solution by backward substitution as follows: where i = 1, 2,, n. ˆ Time complexity: O(n 3 ). x i = b i n j=i+1 ā ij x j, Zheng-Liang Lu 329 / 343

34 Exercise 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 for i = 1 : 3 8 for j = i : 3 9 b(j) = b(j) / A(j, i); % why first? 10 A(j, :) = A(j, :) / A(j, i); 11 end 12 for j = i + 1 : 3 13 A(j, :) = A(j, :) - A(i, :); 14 b(j) = b(j) - b(i); 15 end 16 end 17 x = zeros(3, 1); Zheng-Liang Lu 330 / 343

35 18 for i = 3 : -1 : 1 19 x(i) = b(i); 20 for j = i + 1 : 1 : 3 21 x(i) = x(i) - A(i, j) * x(j); 22 end 23 end 24 x Zheng-Liang Lu 331 / 343

36 More Functions of Linear Algebra 19 ˆ Matrix properties: norm, rcond, det, null, orth, rank, rref, trace, subspace. ˆ Matrix decomposition: lu, chol, qr. 19 See Zheng-Liang Lu 332 / 343

37 Example: 2D Laplace s Equation ˆ A partial differential equation (PDE) is a differential equation that contains unknown multivariable functions and their partial derivatives. 20 ˆ Let Φ(x, y) be a scalar field on R 2. ˆ Consider Laplace s equation 21 as follows: 2 Φ(x, y) = 0, where 2 = is the Laplace operator. x 2 y 2 ˆ Consider the system shown in the next page. 20 See 21 Pierre-Simon Laplace ( ). Zheng-Liang Lu 333 / 343

38 V5 1 V10 V15 V20 V V4 V9 V14 V19 V V3 0.5 V8 V13 V18 V V2 V7 V12 V17 V V1 V6 V11 V16 V ˆ 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 = 1. Zheng-Liang Lu 334 / 343

39 An Simple Approximation 22 ˆ As you can see, we partition the region into many subregions by applying a proper mesh generation. ˆ Then Φ(x, y) can be approximated by Φ(x, y) Φ(x + h, y) + Φ(x h, y) + Φ(x, y + h) + Φ(x, y h), 4 where h is small enough. 22 See _The_Laplace_operator. Zheng-Liang Lu 335 / 343

40 Matrix Formation ˆ By collecting all constraints, we have Ax = b where A = and b = [ ] T. Zheng-Liang Lu 336 / 343

41 Dimension Reduction by Symmetry ˆ As you can see, V 7 = V 17, V 8 = V 18 and V 9 = V 19. ˆ So we can reduce A to A and A = b = [ ] T. ˆ The dimensions of this problem are cut to 6 from 9. ˆ This trick helps to alleviate the curse of dimensionality See Zheng-Liang Lu 337 / 343

42 V V V V V V4 V9 V14 V19 V V V V V V V2 V7 V12 V17 V V1 V6 V11 V16 V Zheng-Liang Lu 338 / 343

43 Remarks ˆ This is a toy example for numerical methods of PDEs. ˆ We can use the PDE toolbox for this case. (Try.) ˆ You may consider the finite element method (FEM). 24 ˆ The mesh generation is also crucial for numerical methods. 25 ˆ You can use the toolbox of Computational Geometry for the triangle mesh See 25 See 26 See https: // Zheng-Liang Lu 339 / 343

44 Example: Method of Least Squares 27 ˆ The method of least squares is a standard approach to the approximate solution of overdetermined systems (m > n). ˆ 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 residuals estimates associated with the data is given by n S = ɛ 2 i. ˆ The best fit in the least-squares sense minimizes the sum of squared residuals. i=1 27 See Zheng-Liang Lu 340 / 343

45 Zheng-Liang Lu 341 / 343

46 Linear Least Squares ˆ The approach is called linear least squares since the assumed function is linear in the parameters to be estimated. ˆ For example, consider y = ax + b, where a and b are to be determined. ˆ Then we have ɛ i = (ax i + b) ŷ i so that S = n ((ax i + b) ŷ i ) 2. i=1 Zheng-Liang Lu 342 / 343

47 ˆ Now consider the partial derivatives of S with respective to a and b: 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 ˆ We reorganize the above equations as follows: n n n a + b x i = x i y i, i=1 x 2 i a i=1 n x i + nb = i=1 i=1 n y i. i=1 ˆ It could be represented by normal equations See (mathematics)#derivation_of_the_normal_equations. Zheng-Liang Lu 343 / 343

Introduction to Matlab Programming with Applications

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