草庐IT

【最小二乘法】矩阵法解法及Python代码

Lil_Yau 2023-10-25 原文

当线性方程组无解时,如何寻求一个最精确解,以贴合在坐标系中离散的坐标点。
给出三个笛卡尔坐标系下的坐标点:
( 1 , 0 )( 2 , 0 )( 1 , 2 ) (1,0)(2,0)(1,2) 1,0)(2,0)(1,2
找到一条直线 y = k x + b y=kx+b y=kx+b,能够贴合这三点(显然不可能)。为表述方便,以 x 1 x₁ x1 k k k x 2 x₂ x2 b b b
得到以下方程组:
{ 1 × k + b = 0 2 × k + b = 0 1 × k + b = 2 ⇔ { 1 × x 1 + x 2 = 0 2 × x 1 + x 2 = 0 1 × x 1 + x 2 = 2 \left\{ \begin{array}{c} 1×k+b=0 \\ 2×k+b=0 \\ 1×k+b=2 \end{array} \right. ⇔ \left\{ \begin{array}{c} 1×x_1+x_2=0 \\ 2×x_1+x_2=0 \\ 1×x_1+x_2=2 \end{array} \right. 1×k+b=02×k+b=01×k+b=2 1×x1+x2=02×x1+x2=01×x1+x2=2
将线性方程组写成矩阵形式:
A = [ 1 1 2 1 1 1 ]   b = [ 0 0 2 ] \begin{gathered} A= \begin{bmatrix} 1&1 \\ 2&1 \\ 1&1 \end{bmatrix} \end{gathered} \begin{gathered}   b= \begin{bmatrix} 0 \\ 0 \\ 2 \end{bmatrix} \end{gathered} A= 121111 b= 002
得出矩阵符合关系式: A x = b Ax=b Ax=b,而我们的目的则是求得一个原始列向量:
x = [ x 1 x 2 ] \begin{gathered} x= \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \end{gathered} x=[x1x2]

分析线性方程组解的情况

(1)

根据矩阵 A = [ 1 1 2 1 1 1 ] \begin{gathered} A= \begin{bmatrix} 1&1 \\ 2&1 \\ 1&1 \end{bmatrix} \end{gathered} A= 121111 列向量线性无关,可得 m = 3 , n = 2 , r = 2 m=3,n=2,r=2 m=3,n=2,r=2
r < m , 且 n = r r<m,且n=r r<m,n=r可知:线性方程组无解,或存在唯一解。

(2)

由于向量b= [ 0 0 2 ] \begin{gathered} \begin{bmatrix} 0 \\ 0 \\ 2 \end{bmatrix} \end{gathered} 002 ,无法由 A = [ 1 1 2 1 1 1 ] \begin{gathered} A= \begin{bmatrix} 1&1 \\ 2&1 \\ 1&1 \end{bmatrix} \end{gathered} A= 121111 线性组合得到,即向量b不在列空间C(A)上。

综上所述可知,线性方程组无解,准确来说,无精确解。

寻找近似解

由线性方程组 A x = b ,即可以得到矩阵 Ax=b,即可以得到矩阵 Ax=b,即可以得到矩阵 A = [ 1 1 2 1 1 1 ] \begin{gathered} A= \begin{bmatrix} 1&1 \\ 2&1 \\ 1&1 \end{bmatrix} \end{gathered} A= 121111 的两个线性无关向量:
a 1 = [ 1 2 1 ] a 2 = [ 1 1 1 ] \begin{gathered} a_1= \begin{bmatrix} 1 \\ 2\\1 \end{bmatrix} \end{gathered} \begin{gathered} a_2= \begin{bmatrix} 1 \\ 1\\1 \end{bmatrix} \end{gathered} a1= 121 a2= 111
这两个向量张成 R 3 R^3 R3空间中的一个二维空间,即矩阵A的列空间 C ( A ) C(A) CA。由于向量 b = [ 0 0 2 ] \begin{gathered} b= \begin{bmatrix} 0 \\ 0 \\ 2 \end{bmatrix} \end{gathered} b= 002 在列空间 C ( A ) C(A) CA外,因此线性方程组 A x = b Ax=b Ax=b是无解的。

虽然线性方程组 A x = b Ax=b Ax=b不存在精确解,但我们可以寻求距离目标最近的近似解。根据上图,我们可以直观的思考到一个合理的方法——从向量 b b b向二维空间即列空间 C ( A ) C(A) C(A)上引垂线,得到 b b b的投影向量 p p p。由 p p p b b b可得误差向量 e = b − p e=b-p e=bp
向量 b b b经“某一矩阵”变换后,可得位于列空间 C ( A ) C(A) C(A)上的投影向量 p p p。而“某一矩阵”则称为“投影矩阵”。我们可以得到关系式:
p = P b p=Pb p=Pb
引出投影向量 p p p后,问题就清晰了,我们可以求出 p p p向量来近似求解方程 A x ^ = p A\hat{x}=p Ax^=p对应 x = [ x 1 x 2 ] \begin{gathered} x= \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \end{gathered} x=[x1x2]

求解的思路与步骤

由上图可知 p p p向量的两个性质:
(1)
p p p向量位于列空间 C ( A ) C(A) CA上,即 A x ^ = p A\hat{x}=p Ax^=p
A x ^ = [ a 1 a 2 ] [ x 1 ^ x 2 ^ ] = p A\hat{x}=\begin{gathered} \begin{bmatrix} a_1&a_2\end{bmatrix} \end{gathered}\begin{gathered} \begin{bmatrix} \hat{x_1} \\ \hat{x_2} \end{bmatrix} \end{gathered}=p Ax^=[a1a2][x1^x2^]=p
(2)
误差向量垂直于二维空间 C ( A ) C(A) C(A)的任一组向量。
a 1 ・ e = 0 和 a 2 ・ e = 0 a_1・e=0和a_2・e=0 a1e=0a2e=0
{ e = b − p p = A x ^ \left\{ \begin{array}{c} e=b-p \\ p=A\hat{x} \end{array} \right. {e=bpp=Ax^代入上述等式:
{ a 1 e = a 1 ( b − p ) = a 1 ( b − A x ^ ) = 0 a 2 e = a 2 ( b − p ) = a 2 ( b − A x ^ ) = 0 \left \{ \begin{array}{c} a_1e=a_1(b-p)=a_1(b-A\hat{x})=0 \\ a_2e=a_2(b-p)=a_2(b-A\hat{x})=0 \end{array} \right. {a1e=a1(bp)=a1(bAx^)=0a2e=a2(bp)=a2(bAx^)=0

将向量内积转化为矩阵乘法

{ a 1 T ( b − A x ^ ) = 0 a 2 T ( b − A x ^ ) = 0 \left \{ \begin{array}{c} a^T_1(b-A\hat{x})=0 \\ a^T_2(b-A\hat{x})=0 \end{array} \right. {a1T(bAx^)=0a2T(bAx^)=0
整理得:
[ a 1 T a 2 T ] ( b − A x ^ ) = 0 \begin{gathered} \begin{bmatrix} a^T_1 \\ a^T_2 \end{bmatrix}(b-A\hat{x})=0 \end{gathered} [a1Ta2T](bAx^)=0
A = [ a 1 a 2 ] 则 [ a 1 T a 2 T ] = A T ⇒ \begin{gathered} A=\begin{bmatrix} a_1 & a_2 \end{bmatrix} \end{gathered}则\begin{gathered} \begin{bmatrix} a^T_1 \\ a^T_2 \end{bmatrix}=A^T \end{gathered}⇒ A=[a1a2][a1Ta2T]=AT
x ^ = ( A T A ) − 1 A T b \hat{x}=(A^TA)^{-1}A^Tb x^=(ATA)1ATb
将上式代入 p = A x ^ p=A\hat{x} p=Ax^,得:
p = A ・ ( A T A ) − 1 A T b p=A・(A^TA)^{-1}A^Tb p=A(ATA)1ATb
代入 p = P b p=Pb p=Pb,得:
P = A ・ ( A T A ) − 1 A P=A・(A^TA)^{-1}A P=A(ATA)1A
公式汇总:
{ x ^ = ( A T A ) − 1 A T b p = A ・ ( A T A ) − 1 A T b P = A ・ ( A T A ) − 1 A \left\{ \begin{array}{c} \hat{x}=(A^TA)^{-1}A^Tb \\ p=A・(A^TA)^{-1}A^Tb \\ P=A・(A^TA)^{-1}A \end{array} \right. x^=(ATA)1ATbp=A(ATA)1ATbP=A(ATA)1A

结论

最小二乘法的几何意义是高维空间中的一个向量在低维空间上的投影。

Python代码求解

代码

import numpy as np
from scipy import linalg
A=np.array([[1,1],
            [2,1],
            [1,1]])
b=np.array([0,0,2])
A_T_A=np.dot(A.T,A)
A_T_A_in=linalg.inv(A_T_A)
x_hat=np.dot(np.dot(A_T_A_in,A.T),b)
p=np.dot(A,x_hat)
p_matrix=np.dot(A,np.dot(A_T_A_in,A.T))
print("近似解x为:")
print(x_hat)
print("\n投影向量p为:")
print(p)
print("\n投影矩阵P为:")
print(p_matrix)

运行结果

近似解x为:
[-1.  2.]    

投影向量p为:
[1. 0. 1.]   

投影矩阵P为:
[[0.5 0.  0.5]
 [0.  1.  0. ]
 [0.5 0.  0.5]]

有关【最小二乘法】矩阵法解法及Python代码的更多相关文章

  1. python - 如何使用 Ruby 或 Python 创建一系列高音调和低音调的蜂鸣声? - 2

    关闭。这个问题是opinion-based.它目前不接受答案。想要改进这个问题?更新问题,以便editingthispost可以用事实和引用来回答它.关闭4年前。Improvethisquestion我想在固定时间创建一系列低音和高音调的哔哔声。例如:在150毫秒时发出高音调的蜂鸣声在151毫秒时发出低音调的蜂鸣声200毫秒时发出低音调的蜂鸣声250毫秒的高音调蜂鸣声有没有办法在Ruby或Python中做到这一点?我真的不在乎输出编码是什么(.wav、.mp3、.ogg等等),但我确实想创建一个输出文件。

  2. ruby - 如何在 buildr 项目中使用 Ruby 代码? - 2

    如何在buildr项目中使用Ruby?我在很多不同的项目中使用过Ruby、JRuby、Java和Clojure。我目前正在使用我的标准Ruby开发一个模拟应用程序,我想尝试使用Clojure后端(我确实喜欢功能代码)以及JRubygui和测试套件。我还可以看到在未来的不同项目中使用Scala作为后端。我想我要为我的项目尝试一下buildr(http://buildr.apache.org/),但我注意到buildr似乎没有设置为在项目中使用JRuby代码本身!这看起来有点傻,因为该工具旨在统一通用的JVM语言并且是在ruby中构建的。除了将输出的jar包含在一个独特的、仅限ruby​​

  3. ruby-on-rails - Rails 源代码 : initialize hash in a weird way? - 2

    在rails源中:https://github.com/rails/rails/blob/master/activesupport/lib/active_support/lazy_load_hooks.rb可以看到以下内容@load_hooks=Hash.new{|h,k|h[k]=[]}在IRB中,它只是初始化一个空哈希。和做有什么区别@load_hooks=Hash.new 最佳答案 查看rubydocumentationforHashnew→new_hashclicktotogglesourcenew(obj)→new_has

  4. ruby-on-rails - 浏览 Ruby 源代码 - 2

    我的主要目标是能够完全理解我正在使用的库/gem。我尝试在Github上从头到尾阅读源代码,但这真的很难。我认为更有趣、更温和的踏脚石就是在使用时阅读每个库/gem方法的源代码。例如,我想知道RubyonRails中的redirect_to方法是如何工作的:如何查找redirect_to方法的源代码?我知道在pry中我可以执行类似show-methodmethod的操作,但我如何才能对Rails框架中的方法执行此操作?您对我如何更好地理解Gem及其API有什么建议吗?仅仅阅读源代码似乎真的很难,尤其是对于框架。谢谢! 最佳答案 Ru

  5. ruby - 模块嵌套代码风格偏好 - 2

    我的假设是moduleAmoduleBendend和moduleA::Bend是一样的。我能够从thisblog找到解决方案,thisSOthread和andthisSOthread.为什么以及什么时候应该更喜欢紧凑语法A::B而不是另一个,因为它显然有一个缺点?我有一种直觉,它可能与性能有关,因为在更多命名空间中查找常量需要更多计算。但是我无法通过对普通类进行基准测试来验证这一点。 最佳答案 这两种写作方法经常被混淆。首先要说的是,据我所知,没有可衡量的性能差异。(在下面的书面示例中不断查找)最明显的区别,可能也是最著名的,是你的

  6. ruby - 寻找通过阅读代码确定编程语言的ruby gem? - 2

    几个月前,我读了一篇关于ruby​​gem的博客文章,它可以通过阅读代码本身来确定编程语言。对于我的生活,我不记得博客或gem的名称。谷歌搜索“ruby编程语言猜测”及其变体也无济于事。有人碰巧知道相关gem的名称吗? 最佳答案 是这个吗:http://github.com/chrislo/sourceclassifier/tree/master 关于ruby-寻找通过阅读代码确定编程语言的rubygem?,我们在StackOverflow上找到一个类似的问题:

  7. ruby - Net::HTTP 获取源代码和状态 - 2

    我目前正在使用以下方法获取页面的源代码:Net::HTTP.get(URI.parse(page.url))我还想获取HTTP状态,而无需发出第二个请求。有没有办法用另一种方法做到这一点?我一直在查看文档,但似乎找不到我要找的东西。 最佳答案 在我看来,除非您需要一些真正的低级访问或控制,否则最好使用Ruby的内置Open::URI模块:require'open-uri'io=open('http://www.example.org/')#=>#body=io.read[0,50]#=>"["200","OK"]io.base_ur

  8. Python 相当于 Perl/Ruby ||= - 2

    这个问题在这里已经有了答案:关闭10年前。PossibleDuplicate:Pythonconditionalassignmentoperator对于这样一个简单的问题表示歉意,但是谷歌搜索||=并不是很有帮助;)Python中是否有与Ruby和Perl中的||=语句等效的语句?例如:foo="hey"foo||="what"#assignfooifit'sundefined#fooisstill"hey"bar||="yeah"#baris"yeah"另外,类似这样的东西的通用术语是什么?条件分配是我的第一个猜测,但Wikipediapage跟我想的不太一样。

  9. java - 什么相当于 ruby​​ 的 rack 或 python 的 Java wsgi? - 2

    什么是ruby​​的rack或python的Java的wsgi?还有一个路由库。 最佳答案 来自Python标准PEP333:Bycontrast,althoughJavahasjustasmanywebapplicationframeworksavailable,Java's"servlet"APImakesitpossibleforapplicationswrittenwithanyJavawebapplicationframeworktoruninanywebserverthatsupportstheservletAPI.ht

  10. 程序员如何提高代码能力? - 2

    前言作为一名程序员,自己的本质工作就是做程序开发,那么程序开发的时候最直接的体现就是代码,检验一个程序员技术水平的一个核心环节就是开发时候的代码能力。众所周知,程序开发的水平提升是一个循序渐进的过程,每一位程序员都是从“菜鸟”变成“大神”的,所以程序员在程序开发过程中的代码能力也是根据平时开发中的业务实践来积累和提升的。提高代码能力核心要素程序员要想提高自身代码能力,尤其是新晋程序员的代码能力有很大的提升空间的时候,需要针对性的去提高自己的代码能力。提高代码能力其实有几个比较关键的点,只要把握住这些方面,就能很好的、快速的提高自己的一部分代码能力。1、多去阅读开源项目,如有机会可以亲自参与开源

随机推荐