JuliaHomotopyContinuation / HomotopyContinuation.jl

A Julia package for solving systems of polynomials via homotopy continuation.
https://www.JuliaHomotopyContinuation.org
MIT License
181 stars 30 forks source link

LinearSubspace produces an unexpected output #461

Closed matthiashimmelmann closed 3 years ago

matthiashimmelmann commented 3 years ago

It seems that there is a problem with the method LinearSubspace. The example from the documentation produces an odd output at best. I double-checked with @alexheaton2 that it is not again an issue with my operating system, as he produces the same output.

julia> A = LinearSubspace([1 0 3; 2 1 3], [5, -2])
1-dim. affine linear subspace {x | Ax=b} with eltype Float64:
A:
[-0.45201106246651834 -0.1594633575218069 -0.8776431148341343; 0.5674520146070269 0.7077338674117247 -0.4208455584140932]
b:
[-0.3464950220709086, -5.65554633590556]
PBrdng commented 3 years ago

Hi!

this is not a real issue. In fact, the documentation is not up to date. (we should fix it!)

The LinearSubspace function automatically orthonormalizes the rows of the matrix; see the documentation.

We can double check your output:

using HomotopyContinuation, LinearAlgebra
A = [1 0 3; 2 1 3.0]
b = [5, -2.0]
L = LinearSubspace(A, b)

A_L = L.extrinsic.A
b_L = L.extrinsic.b
# now A_L and b_L are the matrix/vector from your output.

x₀ = A\b
# x₀ satisfies Ax₀=b

Now:

julia> norm(A_L * x₀ - b_L)
8.899114524108741e-16

and

julia> rank([nullspace(A) nullspace(A_L)])
1

This proves that A x=b and A_L x=b_L describe the same linear spaces.

matthiashimmelmann commented 3 years ago

Thanks, now it makes sense to me! I guess I overread the orthonormality.