14

For a given vector I would like to find the orthogonal basis around it, i.e. the given vector normalized and randomly chosen basis of orthogonal sub-space. Is there a convenient function for this in Julia?

xealits
  • 4,224
  • 4
  • 27
  • 36

3 Answers3

9

The function you are looking for is called nullspace.

julia> x = randn(5);

julia> x⊥ = nullspace(x');

julia> x'x⊥
1×4 Array{Float64,2}:
 7.69373e-16  -5.45785e-16  -4.27252e-17  1.26778e-16
Andreas Noack
  • 1,370
  • 6
  • 11
  • Indeed! (Also, pretty obvious..) And with the normalized `x` the full basis will be: `hcat(x / sqrt(x'x), nullspace(x'))'` – xealits Sep 20 '16 at 14:17
  • 2
    `x'` used to create a row matrix and it now just creates a `Adjoint` of the vector. You can make it work on recent versions of Julia with `nullspace(Matrix(x'))`. An attempt to make `nullspace(x')` work can be found in https://github.com/JuliaLang/julia/pull/33385/ – Benoît Legat Sep 27 '19 at 10:14
8

You could define a function orth (if someonehasn't already done this)

orth(M) = qr(M)[1]

See here: https://groups.google.com/forum/#!topic/julia-users/eG6a4tj7LGg and http://docs.julialang.org/en/release-0.4/stdlib/linalg/

Or from IterativeSolvers.jl:

orthogonalize{T}(v::Vector{T}, K::KrylovSubspace{T})

See: https://github.com/JuliaMath/IterativeSolvers.jl

Alexander Morley
  • 4,111
  • 12
  • 26
  • 1
    I seriously doubt you're going back to medicine after the PhD, Alex :D (where excel sheets are considered IT wizardry) – Tasos Papastylianou Sep 19 '16 at 23:40
  • 1
    thanks for suggestions! but I'll accept the answer with `nullspace` -- it's pretty basic, but fits the problem perfectly and should not be forgotten – xealits Sep 20 '16 at 14:27
2

The following will calculate an orthogonal basis for matrix M

function orth(M::Matrix)
  matrixRank = rank(M)
  Ufactor = svdfact(M)[:U]
  return Ufactor[:,1:matrixRank]
end

With julia documentation:

"""
orth(M)

Compute an orthogonal basis for matrix `A`.

Returns a matrix whose columns are the orthogonal vectors that constitute a basis for the range of A.
If the matrix is square/invertible, returns the `U` factor of `svdfact(A)`, otherwise the first *r* columns of U, where *r* is the rank of the matrix.
# Examples
```julia
julia> orth([1 8 12; 5 0 7])
2×2 Array{Float64,2}:
 -0.895625  -0.44481
 -0.44481    0.895625
```
```
julia> orth([1 8 12; 5 0 7 ; 6 4 1])
3×3 Array{Float64,2}:
 -0.856421   0.468442   0.217036
 -0.439069  -0.439714  -0.783498
 -0.27159   -0.766298   0.582259
```
"""
function orth(M::Matrix)
  matrixRank = rank(M)
  Ufactor = svdfact(M)[:U]
  return Ufactor[:,1:matrixRank]
end
RaidenF
  • 3,411
  • 4
  • 26
  • 42