在Julia中,两个有理值矩阵的矩阵分割返回一个浮点矩阵。我怎样才能得到一个有理数的矩阵呢?Julia中的有理矩阵分割
我不能只使用convert(Array{Rational}, A \ b)
,因为与浮点数相关的精度损失。
在Julia中,两个有理值矩阵的矩阵分割返回一个浮点矩阵。我怎样才能得到一个有理数的矩阵呢?Julia中的有理矩阵分割
我不能只使用convert(Array{Rational}, A \ b)
,因为与浮点数相关的精度损失。
你必须实现一个枢轴QR factorization算法的有理矩阵,这是一个非常平凡的任务,虽然当然不是不可能的。 Julia使用LAPACK DGEQP3 function为64位浮点矩阵执行此操作。即使你设法得到一个合理的QR分解工作,它也不会像LAPACK算法那么快。所以我想我会问你需要什么精确的有理矩阵分割。
另外:你可能会发现在julia-users list上问这样的问题更有成效,因为你可以谈论它 - 这不是真正的“问题和回答”类型的问题。
更新:现在这“只是工作”,因为通用摆动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
需要注意的是四元数乘法是不可交换的,这使得这相当有趣。
这只是一个感兴趣的问题。我想用Julia来推导高阶积分方法,其中系数是线性系统“Ax = b”的解。精确的有理系数比浮点型更适合数学分析。 –
拥有BLAS/LAPACK功能的相当大的子集的合理类似物当然是一个有趣的项目,但我不知道任何人都承认过。我不知道任何现有的系统可以做这种事情,除了Mathematica - 虽然我没有尝试过,所以我不能证实它。理性溢出可能是这类事业的主要问题。 – StefanKarpinski
我已将此问题发布到此处的julia用户:https://groups.google.com/forum/#!topic/julia-users/RLqRiHm-MpA –