在 Julia 上使用高斯乔丹消元法求矩阵逆

Question

我想制作一种使用高斯-乔丹方法来计算矩阵逆的算法。我对Matlab比较熟悉，所以先写了这个Matlab算法。它工作得很好，但我无法让它在 Julia 上工作

function [Inv,P]=GJinv(A)
  format rat
  n=length(A);
  P=eye(n); 
  Inv=[A,eye(n)]; 
  for k=1:n
        
        [~,r] = max(abs(A(k:end,k)));
        r = n-(n-k+1)+r;    
        #αλλάζω γραμμές#
        Inv([k r],:) = Inv([r k],:);
        P([k r],:) = P([r k],:);
       
        Inv(k,:) = Inv(k,:) / Inv(k,k);
        for i=1:n
          if i~=k
            Inv(i,:)=Inv(i,:)-(Inv(i,k)*Inv(k,:)); 
          endif
        endfor
  endfor
  Inv=Inv(1:n,n+1:2*n);
endfunction

到目前为止我的工作：

using LinearAlgebra
function GJinv(A)
    (n,b)=size(A)
    Inv=[A,I(n)] 
    for k=1:n
          
          local ma=maximum(abs.(A[k:n,k:k])) 
          local pos = findfirst( x -> x == ma, abs.(A[k:n,k:n]))
          local r=pos[1] ; local r= n-(n-k+1)+r 
          #αλλάζω γραμμές#
          Inv[[k, r],:] = Inv[[r, k],:]
          
          Inv[k,:] = Inv[k,:] / Inv[k,k]
          for i=1:n
            if i!=k
              Inv[i,:]=Inv[i,:]-(Inv[i,k]*Inv[k,:]) 
            end
          end
    end
    Inv=Inv[1:n,n+1:2*n]
    print(Inv)
end
#TEST
A=[2 2 3; 4 5 6; 1 2 4]
GJinv(A)    

我收到以下错误：

ERROR: LoadError: InexactError: Int64(1.25)
Stacktrace:
  [1] Int64
    @ .\float.jl:912 [inlined]
  [2] convert
    @ .\number.jl:7 [inlined]
  [3] setindex!
    @ .\array.jl:1024 [inlined]
  [4] macro expansion
    @ .\multidimensional.jl:960 [inlined]
  [5] macro expansion
    @ .\cartesian.jl:64 [inlined]
  [6] _unsafe_setindex!(::IndexLinear, ::Matrix{Int64}, ::Vector{Float64}, ::Int64, ::Base.Slice{Base.OneTo{Int64}})
    @ Base .\multidimensional.jl:955
  [7] _setindex!
    @ .\multidimensional.jl:944 [inlined]
  [8] setindex!
    @ .\abstractarray.jl:1396 [inlined]
  [9] GJinv(A::Matrix{Int64})
    @\GJinv.jl:13
 [10] top-level scope
    @\GJinv.jl:24
in expression starting at @\GJinv.jl:24

任何帮助将不胜感激。

Answer 1

你被两个问题所困扰：MATLAB 的语法与 Julia 的语法不太一样。

首先，如果您在 Julia 中声明一个包含所有整数值的矩阵，则该矩阵将具有 Int 类型，这将在您稍后除法时给出问题，因为在 Julia 中 /(Int x, Int y) 的类型是浮点 x / 是。

其次，正如 S Zakharov 上面提到的，Inv[A, I(n)] 会给你一个 Julia 中的矩阵向量，但你想要由两个矩阵组成的水平 (hcat) 矩阵。

所以考虑这个修正（我修正了语法，但没有改变你的代数afaik）：

使用线性代数

function GJinv(squaremat)
    n, b = size(squaremat)
    @assert n == b # number of rows must == number of cols
    A = float.(squaremat)
    Inv = hcat(A, I(n))
    for k = 1:n
        ma = maximum(abs.(A[k:n,k:k])) 
        r = k + findfirst( x -> x == ma, abs.(A[k:n,k:n]))[1] - 1
        #αλλάζω γραμμές#
        Inv[[k, r],:] = Inv[[r, k],:]
        Inv[k,:] = Inv[k,:] / Inv[k,k]
        for i=1:n
            if i!=k
                Inv[i,:]=Inv[i,:]-(Inv[i,k]*Inv[k,:]) 
            end
        end
    end
    Inv = round.(Inv[1:n,n+1:2*n], digits = 6) # rounding makes it print better
    print(Inv)
end

 #TEST
 A=[2 2 3; 4 5 6; 1 2 4]

 GJinv(A)  # [1.6 -0.4 -0.6; -2.0 1.0 0.0; 0.6 -0.4 0.4]