【问题标题】:Solving underdetermined linear systems with R用 R 求解欠定线性系统
【发布时间】:2018-10-26 13:08:33
【问题描述】:

R 可以解决欠定线性系统:

A = matrix((1:12)^2,3,4,T)
B = 1:3
qr(A)$rank  # 3
qr.solve(A, B)  # solutions will have one zero, not necessarily the same one
# 0.1875 -0.5000  0.3125  0.0000
solve(qr(A, LAPACK = TRUE), B)
# 0.08333333 -0.18750000  0.00000000  0.10416667

(它给出了无数解中的一种解)。

但是,如果排名(此处为 2)低于行数(此处为 3),则不起作用:

A = matrix(c((1:8)^2,0,0,0,0),3,4,T)
B = c(1,2,0)
A
#      [,1] [,2] [,3] [,4]
# [1,]    1    4    9   16
# [2,]   25   36   49   64
# [3,]    0    0    0    0

qr.solve(A, B)  # Error in qr.solve(A, B) : singular matrix
solve(qr(A, LAPACK = TRUE), B)  # Error in qr.coef(a, b) : error code 3 

但是这个系统确实有解决方案!

我知道一般的解决方案是使用 SVD 或 A 的广义/伪逆(参见this question 及其答案),但是:

solveqr.solve 是否存在自动将系统 AX=B 减少到仅具有 rank(A) 行的等效系统 CX=D 的方法,其中qr.solve(C, D)开箱即用?

例子:

C = matrix(c((1:8)^2),2,4,T)
D = c(1,2)
qr.solve(C, D)
# -0.437500  0.359375  0.000000  0.000000

【问题讨论】:

    标签: r math linear-algebra numerical-methods


    【解决方案1】:

    qr.coefqr 似乎可以完成这项工作:

    (A <- matrix(c((1:8)^2, 0, 0, 0, 0), nrow = 3, ncol = 4, byrow = TRUE))
    #     [,1] [,2] [,3] [,4]
    # [1,]    1    4    9   16
    # [2,]   25   36   49   64
    # [3,]    0    0    0    0
    (B <- c(1, 2, 0))
    # [1] 1 2 0
    (X0 <- qr.coef(qr(A), B))
    # [1] -0.437500  0.359375        NA        NA
    X0[is.na(X0)] <- 0
    X0
    # [1] -0.437500  0.359375  0.000000  0.000000
    # Verification:
    A %*% X0
    #      [,1]
    # [1,]    1
    # [2,]    2
    # [3,]    0
    

    第二个例子:

    (A<-matrix(c(1, 2, 0, 0, 1, 2, 0, 0, 1, 2, 1, 0), nrow = 3, ncol = 4, byrow = TRUE))
    #      [,1] [,2] [,3] [,4]
    # [1,]    1    2    0    0
    # [2,]    1    2    0    0
    # [3,]    1    2    1    0
    (B<-c(1, 1, 2))
    # [1] 1 1 2
    qr.solve(A, B)
    # Error in qr.solve(A, B) : singular matrix 'a' in solve
    (X0 <- qr.coef(qr(A), B))
    # [1]  1 NA  1 NA
    X0[is.na(X0)] <- 0
    X0
    # [1] 1 0 1 0
    A %*% X0
    #      [,1]
    # [1,]    1
    # [2,]    1
    # [3,]    2
    

    【讨论】:

    • 是的,但在这里您必须手动设置要保留的行 (1:2)。如果系统有 100 行并且排名为 70,即必须丢弃 30(可能不连续)行怎么办? (或者我在你的例子中遗漏了一些东西)。
    • 重新表述我的问题@JuliusVainora:如何创建一个解决方案 X0 使得 A X0 = B,而不指定应该保留哪些行以及应该丢弃哪些行? (代码中的“1:2”)
    • 没太注意,但是w &lt;- !is.na(coef); A[,w] %*% na.omit(coef)有用吗?
    • 没有必要知道qr.coef(qr(A), B) 中的1:2,这就是您的全部问题。我只添加了一个验证步骤。至于摆脱NA,您可以使用coef[!is.na(coef)],或者@BenBolker 建议的验证部分。
    • @JuliusVainora qr.coef(qr(A), B) 究竟代表什么?问题是关于在这个例子中获得一个只有两行的等效系统CX=D。如何从qr.coef(qr(A), B) 到矩阵CD?我看到我们已经很接近了,但我只是想确保使用最好的方法。
    猜你喜欢
    • 1970-01-01
    • 2013-11-14
    • 2017-12-13
    • 2018-01-17
    • 2020-04-03
    • 2019-01-15
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多