matlab - JULIA 中无需旋转的 LU 分解

Question

尝试根据此答案重写 lu_nopivot https://stackoverflow.com/a/41151228进入 JULIA 并仅使用一个循环。 我写的 Julia 代码

using LinearAlgebra
function lu_nopivot(A)
    n = size(A, 1)
    L = Matrix{eltype(A)}(I, n, n)
    U = copy(A)
    for k = 1:n
        L[k+1:n,k] = U[k+1:n,k] / U[k,k]
        U[k+1:n,:] = U[k+1:n,:] - L[k+1:n,k]*U[k,:]  
    end
    return L, U
end

但调用该函数 L, U = lu_nopivot(A) 在 上给出错误 MethodError: no methodmatching *(::Vector{Float64},::Vector{Float64}) L[k+1:n,k]*U[k,:]。这是什么原因？

Answer 1

问题是，当您执行 U[k, :] 时，即使您提取的是一行，您得到的也是一个列向量。因此，L[k+1:n,k]*U[k,:] 变成了将列向量乘以列向量的尝试。

在 Julia 中获取行向量(又名 1 x N 矩阵)的一种方法(尽管我不知道这是否是惯用的方法)是执行 U[k:k , :] 改为:

U[k+1:n,:] = U[k+1:n,:] - L[k+1:n,k] * U[k:k,:]  

但请注意，Julia 的 lu 实现已经有一个 no-pivot 选项:

Pivoting can be turned off by passing pivot = NoPivot().

julia> M = rand(1.0:100.0, 3, 3)
3×3 Matrix{Float64}:
 71.0  50.0  23.0
 82.0  63.0  37.0
 96.0  28.0   5.0

julia> L, U = lu_nopivot(M); # after the above k:k change has been made

julia> L
3×3 Matrix{Float64}:
 1.0       0.0      0.0
 1.15493   1.0      0.0
 1.35211  -7.53887  1.0

julia> U
3×3 Matrix{Float64}:
 71.0  50.0      23.0
  0.0   5.25352  10.4366
  0.0   0.0      52.5818

julia> lu(M, NoPivot())
LU{Float64, Matrix{Float64}}
L factor:
3×3 Matrix{Float64}:
 1.0       0.0      0.0
 1.15493   1.0      0.0
 1.35211  -7.53887  1.0
U factor:
3×3 Matrix{Float64}:
 71.0  50.0      23.0
  0.0   5.25352  10.4366
  0.0   0.0      52.5818

使用它可能会更高效，也更健壮(例如，您当前的实现无法处理 eltype Int 的矩阵，因为 L 和 U 被赋予与输入相同的类型，但需要包含Float)。