朱莉娅:交换出错

Question

我正在使用Julia 0.3.4

我正在尝试使用高斯消元法编写LU分解.所以我必须交换行.这是我的问题:

如果我正在使用a,b = b,a我收到错误,

但如果我正在使用:

function swapRows(row1, row2)
    temp = row1
    row1 = row2
    row2 = temp
end

一切都很好.

我做错了什么或者是个错误？

这是我的源代码:

function lu_t(A::Matrix)
  # input value: (A), where A is a matrix
  # return value: (L,U), where L,U are matrices

  function swapRows(row1, row2)
    temp = row1
    row1 = row2
    row2 = temp
    return null
  end

  if size(A)[1] != size(A)[2]
    throw(DimException())
  end

  n = size(A)[1] # matrix dimension
  U = copy(A) # upper triangular matrix
  L = eye(n) # lower triangular matrix

  for k = 1:n-1 # direct Gaussian elimination for each column `k`
    (val,id) = findmax(U[k:end,k]) # find max pivot element and it's row `id`

    if val == 0 # check matrix for singularity
      throw(SingularException())
    end

    swapRows(U[k,k:end],U[id,k:end]) # swap row `k` and `id`
    # U[k,k:end],U[id,k:end] = U[id,k:end],U[k,k:end] - error

    for i = k+1:n # for each row `i` > `k`
      ? = U[i,k] / U[k,k] # find elimination coefficient `?`
      L[i,k] = ? # save to an appropriate position in lower triangular matrix `L`
      for j = k:n # update each value of the row `i`
        U[i,j] = U[i,j] - ??U[k,j]
      end
    end
  end

  return (L,U)
end

###### main code ######
A = rand(4,4)
@time (L,U) = lu_t(A)
@test_approx_eq(L*U, A)

Answer 1

该swapRows函数是一个无操作,并且没有任何影响 - 它所做的就是交换一些局部变量名.请参阅有关分配和变异之间差异的各种讨论:

• https://groups.google.com/d/msg/julia-users/oSW5hH8vxAo/llAHRvvFVhMJ
• http://julia.readthedocs.org/en/latest/manual/faq/#i-passed-an-argument-x-to-a-function-modified-it-inside-that-function-but-on-the -outside最可变-X-是静止不变的,为什么
• http://julia.readthedocs.org/en/latest/manual/faq/#why-does-xy-allocate-memory-when-x-and-y-are-arrays

常量null并不意味着你的想法 - 在Julia v0.3中,它是一个计算线性变换的零空间的函数; 在Julia v0.4中它仍然意味着这个但已被弃用并重命名为nullspace.朱莉娅的"无趣"价值被称为nothing.

我不确定您注释掉的行交换代码有什么问题,但这种通用方法确实有效:

julia> X = rand(3,4)
3x4 Array{Float64,2}:
 0.149066  0.706264   0.983477  0.203822
 0.478816  0.0901912  0.810107  0.675179
 0.73195   0.756805   0.345936  0.821917

julia> X[1,:], X[2,:] = X[2,:], X[1,:]
(
1x4 Array{Float64,2}:
 0.478816  0.0901912  0.810107  0.675179,

1x4 Array{Float64,2}:
 0.149066  0.706264  0.983477  0.203822)

julia> X
3x4 Array{Float64,2}:
 0.478816  0.0901912  0.810107  0.675179
 0.149066  0.706264   0.983477  0.203822
 0.73195   0.756805   0.345936  0.821917

由于这会创建一对我们无法消除分配的临时数组,因此这不是最有效的方法.如果你想要最有效的代码,循环两行并交换标量值对会更快:

function swapRows!(X, i, j)
    for k = 1:size(X,2)
        X[i,k], X[j,k] = X[j,k], X[i,k]
    end
end

请注意,在Julia中,通常使用尾随命名函数来改变一个或多个参数!.目前,闭包(即内部函数)存在一些性能问题,因此您需要在顶级作用域中定义这样的辅助函数,而不是以您获得它的方式在另一个函数内部定义.

最后,我认为这是一个练习,因为Julia附带了经过仔细调整的通用(即它适用于任意数字类型)LU分解:http://docs.julialang.org/en/release-0.3/stdlib/linalg/#Base.路.

-