julia - Julia 中的有理矩阵除法

Question

在 Julia 中，两个有理值矩阵的矩阵除法返回一个浮点矩阵。我怎样才能获得有理数矩阵呢？

我不能只用 convert(Array{Rational}, A \ b)因为与浮点数相关的精度损失。

Answer 1

你必须实现一个旋转的 QR factorization有理矩阵的算法，这是一项非常重要的工作，尽管肯定不是不可能的。 Julia 使用 LAPACK DGEQP3 function对 64 位浮点矩阵执行此操作。即使您确实设法使合理的 QR 分解工作，它也远不及 LAPACK 算法那么快。所以我想我会问你需要精确的有理矩阵除法。

旁白:您可能会发现在 julia-users list 上提出这样的问题更有成效。因为您将能够就此进行对话 - 这不是真正的“询问和回答”类型的问题。

更新:这现在“有效”，因为 Julia 中现在存在一个通用的旋转 QR:

julia> A = [rand(1:10)//rand(1:10) for i=1:4, j=1:4]
4x4 Array{Rational{Int64},2}:
 5//3  1//2  10//1   8//9
 1//1  3//2   6//7   2//3
 4//1  8//9   6//7  10//3
 7//2  5//2   1//2   5//1

julia> b = collect(1:4)
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> c = A\b
4-element Array{Rational{Int64},1}:
  42055//62497
  61344//62497
  -2954//62497
 -19635//124994

julia> A*c
4-element Array{Rational{Int64},1}:
 1//1
 2//1
 3//1
 4//1

但是，请注意 Rational{Int}相当容易溢出，因此您可能需要使用 Rational{Int128}或 Rational{BigInt}以避免溢出。这个算法是完全通用的，也适用于更奇特的数值类型:

julia> using Quaternions # install with Pkg.add("Quaternions")

julia> A = [Quaternion(rand(1:10), rand(1:10), rand(1:10), rand(1:10)) for i=1:4, j=1:4]
4x4 Array{Quaternions.Quaternion{Int64},2}:
  4 + 3im + 5jm + 8km  9 + 7im + 10jm + 3km    9 + 3im + 1jm + 7km    8 + 4im + 8jm + 5km
  1 + 4im + 2jm + 1km   5 + 4im + 1jm + 4km    1 + 5im + 8jm + 2km    7 + 2im + 5jm + 3km
 10 + 1im + 4jm + 4km   2 + 4im + 9jm + 2km  2 + 10im + 4jm + 10km    2 + 3im + 4jm + 8km
  7 + 4im + 3jm + 7km   8 + 3im + 5jm + 9km    6 + 5im + 1jm + 3km  10 + 10im + 3jm + 1km

julia> b = collect(1:4)
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> c = A\b
4-element Array{Quaternions.Quaternion{Float64},1}:
    0.18112 + 0.019288355350921868im - 0.002638049486498678jm + 0.11047233503816825km
 0.000578119 + 0.08073035854610844im - 0.023278758601757744jm - 0.16761078955242836km
  -0.0501257 - 0.02428891715971441im - 0.11059096793480247jm - 0.059017235311989824km
   0.0394953 - 0.16771397199827004im - 0.019823106776170954jm + 0.05251791790026253km

julia> A*c
4-element Array{Quaternions.Quaternion{Float64},1}:
                                    1.0 + 1.1102230246251565e-16im + 0.0jm + 0.0km
                                     2.0 + 2.220446049250313e-16im + 0.0jm + 0.0km
 3.0 + 4.440892098500626e-16im - 2.220446049250313e-16jm + 8.881784197001252e-16km
 4.0 + 2.220446049250313e-16im - 2.220446049250313e-16jm + 6.661338147750939e-16km

julia> norm(b - A*c)
1.680072297996111e-15

请注意，四元数乘法不是可交换的，这使得这很有趣。