1 / 29

Least Squares

Least Squares. 给定一组测量或观测值. 假设 y 是 n 个函数的线性组合. 设 m×n 矩阵 X 定义为. 所以. β = Xy. 直线：. 多项式：. 有理函数：. 指数：. Log-linear:. Gaussians:. 中国人口预测. p = 1e4*[ ]'; t = (1949:1:2005)'; x = (1940: 0.2 :2019)'; w = 2010; d = 3; s = (t-1979)/30; c = polyfit(s,p,d);

natara
Download Presentation

Least Squares

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. Least Squares

  2. 给定一组测量或观测值 假设y是n个函数的线性组合 设m×n矩阵X定义为 所以 β=X\y

  3. 直线： 多项式： 有理函数： 指数：

  4. Log-linear: Gaussians:

  5. 中国人口预测 p = 1e4*[ ]'; t = (1949:1:2005)'; x = (1940: 0.2 :2019)'; w = 2010; d = 3; s = (t-1979)/30; c = polyfit(s,p,d); s = (x-1979)/30; y = polyval(c,s); s = (w-1979)/30; z = polyval(c,s); plot(t,p,'o',x,y,'b-',w,z,'r*') s = (t-1979)/30; c = polyfit(s,p,3); X = [s.^3, s.^2, s, ones(length(s),1)]; X\p

  6. 1、求逆的过程困难 2、

  7. 写成矩阵形式为： 其中，

  8. Matlab中避免用法方程求解 QR分解： Q正交阵(orthogonal) R上三角阵 求解：

  9. s: -0.9667 -0.6333 -0.3000 0.0333 0.3667 0.7000 y: 55196 66207 82992 98705 114333 126743 X: 0.9344 -0.9667 1.0000 0.4011 -0.6333 1.0000 0.0900 -0.3000 1.0000 0.0011 0.0333 1.0000 0.1344 0.3667 1.0000 0.4900 0.7000 1.0000

  10. u(1) = u(1) + sigma; rho = 1/( sigma *u(1)); X(i,k) = 0; X(k,k) = -sigma; j = k+1:n; v = rho*(u'* X(i,j)); X(i,j) = X(i,j) - u*v; if ~isempty(y) tau = rho*(u'* y(i)); y(i) = y(i) - tau*u; end end end beta = X(1:n,1:n)\y(1:n) polyval(beta,(2010-1979)/30) s = ( (1950:10:2000)'-1979 )/30; X = [s.^2, s, ones(size(s),1)] y = [ ] [m,n] = size(X); for k = 1:min(m-1,n) i = k:m; u = X(i,k); sigma = norm(u); if sigma ~= 0 if u(1) ~= 0, sigma = sign(u(1))*sigma; end u(1) = u(1) + sigma; rho = 1/( sigma *u(1));

  11. X = -1.1403 0.6945 -1.7987 0 -0.3122 0.4589 0 -0.2279 0.8786 0 0.0342 0.9985 0 0.4743 0.8186 0 1.0923 0.3390 y = 1.0e+005 * -1.4311 0.2787 0.7439 0.9860 1.0148 0.7991 X = -1.1403 0.6945 -1.7987 0 1.2525 0.3587 0 0 0.8640 0 0 1.0007 0 0 0.8490 0 0 0.4090 y = 1.0e+005 * -1.4311 0.9033 0.8349 0.9723 0.8255 0.3630 X = -1.1403 0.6945 -1.7987 0 1.2525 0.3587 0 0 -1.6236 0 0 0 0 0 0 0 0 0 y = 1.0e+005 * -1.4311 0.9033 -1.5667 0.0062 0.0058 -0.0318

  12. 广义逆 列满秩 是左逆 不是右逆,但满足

  13. 秩亏损(Rank Deficiency) 接近秩亏损 X = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 y = 16 17 18 19 20 pinv(X)*y ans = -7.5556 0.1111 7.7778 X\y-pinv(X)*y ans = 0.0556 -0.1111 0.0556 X\y ans = -7.5000 0 7.8333

  14. 非线性最小二乘问题 fminsearch fmincon fminunc lsqnonlin lsqcurvefit

  15. 调试 dbstop at 14 in xxx %F12 dbquit dbtype xxx.m 12:15 %type,edit error(‘xxxxx’) pause(n) keyboard; %K>>return

  16. 分块循环约化（Block cyclic reduction） 求解矩形域上的Poisson方程，产生如下的系数矩阵： 用分块循环约化法求 （假设解全为1） 比较分块循环约化与通常求解方法的效率。 参考 G. H. Golub, C. F. Van Loan, Matrix Computations, 3rd ed., Johns Hopkins University Press, Baltimore and London, 1996. 有中文版，p.202

  17. 进一步探讨用分块循环约化法于下面方程的求解进一步探讨用分块循环约化法于下面方程的求解 Solve 1. Block cyclic reduction (BCR) B. L. Buzbee, G. H. Golub, C.W. Nielson, On direct methods for solving Poisson’s equations, SIAM J. Numer. Anal., 7(1970), pp. 627-656. 2. Block orthogonal factorization (BOF) G. Fairweather, I. Gladwell, Algorithms for almost block diagonal linear systems, SIAM Review, 46(2004), pp. 49-58. S. J. Wright, Stable parallel algorithms for two-point boundary value problem, SIAM J. Sci. Statist. Comput., 13(1992), pp. 742-764.

  18. 不完全LU分解（ILU） 参考Yousef Saad的书 《Iterative methods for sparse linear systems》的第10章. 该书第一版可以在Saad的主页上下载。

  19. Cholesky factorization 1. Basic Cholesky factorization 1) gaxpy Cholesky 2) outer-product Cholesky 3) 分块点积的Cholesky分解 4)其他实现方法 比较各种方法的效率，并与Matlab函数chol()比较。 G. H. Golub, C. F. Van Loan, Matrix Computations, 3rd ed., Johns Hopkins University Press, Baltimore and London, 1996.

  20. 2. Incomplete Cholesky (IC) factorization 参考 G.H. Golub, C. F. Van Loan 的《Matrix Computation》. 中文版 p.620 T. A. Manteuffel, An incomplete factorization technique for positive definite linear systems, Mathematics of Computation, 34, 150(1980), 473-497. J. Meijerink, H. A. Van der Vorst, An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix, Math. Comput., 31, 137(1977), pp. 134-155. R. S. Varga, E. B. Saff, V. Mehrmann, Incomplete factorizations of matrices and connections with H-matrices, SIAM J. Numer. Anal., 17(1980), pp. 787-793.

  21. 3. Robust Incomplete Cholesky factorization 1) T. A. Manteuffel, An incomplete factorization technique for positive definite linear systems, Mathematics of Computation, 34, 150(1980), 473-497. M. Ajiz, A. Jennings, A robust incomplete Choleski-conjugate gradient algorithm, Internat. J. Numer. Methods Engrg., 20, 5(1984), pp. 949-966. 2) M. Tismenetsky, A new preconditioning technique for solving large sparse linear systems, Linear Algebra Appl., 154-156(1991), pp. 331-353. 3) I. E. Kaporin, High quality preconditioning of a general symmetric positive definite matrix based on its -docomposition, Numer. Linear Algebra Appl., 5(1998), pp. 483-509.

  22. Let A be thematrix whose entries are zero everywhere except for the primes 2, 3, 5, 7, …, 224737 along the main diagonal and the number 1 in all the positions with. What is the (1,1) entry of ?

  23. The Lorenz Strange Attractor The Lorenz system is given by Using the parameters and and the initial conditions (1)apply the Runge-Kutta method using step size to solve the system for and then plot each of the projection z vs. x, y vs. x, and y vs. z. (2)Perturb the initial conditions slightly to x(0)=-8.02, y(0)=7.98 (and same z(0)), solve the new system with the same method, and plot the original x minus the new x versus t on the whole time interval.

  24. Numerical algorithms for solving the linear two-point boundary value problem (BVP) • Two approach: • multiple shooting technique • finite-difference method(建议用此法) S. J. Wright, Stable parallel algorithms for two-point boundary value problems, SIAM J. Sci. Stat. Comput., 13(1992), pp. 742-764.

  25. 用Newton法解复数方程 在平面上把收敛到不同根的初始值分别标上不同颜色, 得到Newton法收敛域的彩色图片. 函数fzero的实现 Algorithms on QR算法求特征值问题

More Related