JULIA 中没有旋转的 LU 分解

LU decomposition without pivoting in JULIA

尝试将此答案 中的 lu_nopivot 重写为 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)L[k+1:n,k]*U[k,:] 上给出错误 MethodError: no method matching *(::Vector{Float64}, ::Vector{Float64})。这是什么原因?

问题是,当您执行 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 被赋予与输入相同的类型,但需要包含Floats).