Nemocas / Nemo.jl

Julia bindings for various mathematical libraries (including flint2)
http://nemocas.github.io/Nemo.jl/
Other
185 stars 58 forks source link

Generalized Gram Schmidt Orthogonalization (arb, acb) #734

Open ambiso opened 4 years ago

ambiso commented 4 years ago

Hello!

I was trying to use gso today, however, unfortunately, it's only defined for fmpq_mat, and I wasn't able to use it with a matrix of complex values (acb_mat).

gso simply calls this C function, which doesn't look like a sophisticated implementation to me.

I found this generic version of the Gram Schmidt process here (see section 5.4):

function gram_schmidt(a; tol = 1e-10)
    q = []
    for i = 1:length(a)
        qtilde = a[i]
        for j = 1:i-1
            qtilde -= (q[j]'*a[i]) * q[j]
        end
        if norm(qtilde) < tol
            println("Vectors are linearly dependent.")
            return q
        end
        push!(q, qtilde/norm(qtilde))
    end
    return q
end

Would it be possible to use this or something similar in Nemo?

Kind Regards, ambiso

thofma commented 4 years ago

Not sure how easy this is.

Maybe @fredrik-johansson has an idea.

wbhart commented 4 years ago

You can get massive cancellation in gso, which means a total loss of precision in the worst case.

There's some quite complicated code in the Flint LLL for handling this heuristically, including some algorithms for effectively increasing the precision. Unfortunately, it is only implemented for doubles and mpfr's, not for Arbs.

ambiso commented 4 years ago

If precision is a problem, would it be possible to create a type for arbitrary precision complex number matrices? (using fmpq?)

wbhart commented 4 years ago

Precision isn't the issue. Cancellation is.

wbhart commented 4 years ago

But obviously if you do everything over Q instead of floating point, although very slow, it would work. But we have no type for complex numbers with rational parts. There are currently no plans to add such a type.

You are probably better off using Julia directly for this sort of thing, as it surely has a gso function for complex numbers over Rational{BigInt}. You can certainly create such things in Julia:

julia> complex(BigInt(1)//2, BigInt(2)//3)
1//2 + 2//3*im