zhaowenlan1779 / SVDSketch.jl

Julia implementation of the `svdsketch` function from Matlab
http://blog.zhupengfei.com.cn/SVDSketch.jl/
MIT License
1 stars 0 forks source link

Type instability involving `eigen` and `adjoint` #2

Open ChenNingCong opened 1 year ago

ChenNingCong commented 1 year ago

This package suffers from severe type instability from some LinearAlgebra function. Firstly, it seems that eigen is always performed on a symmetric matrix (A * A'), you should use eigen(Hermitian(x)). Secondly, A' (adjoint) returns a different type instead of Matrix{Float64}, you need to use function barrier to isolate the type instability. Thirdly, alpha in function _svdsketch is initialized with 0, an integer. Then it's assigned with a floating number. You should initialize it with a floating number. These errors are found by SimpleTypeChecker, a fortran-style static type checker for Julia.

zhaowenlan1779 commented 1 year ago

Thank you! I fixed the other two problems, but I'm not quite sure which part you are referring to regarding the adjoint issue.

I also tried to run your checker but it keeps telling me that "AST has no matching green node", not sure how to fix that?

ChenNingCong commented 1 year ago

Technically, you don't need to fix the case of adjoint. Julia's compiler silently removes the type instability for you. A' is shorthand for adjoint(A). But adjoint(A) doesn't return a Matrix{Float64}, it returns a wrapper of different type.

julia> [1 2; 4 5]'
2×2 adjoint(::Matrix{Int64}) with eltype Int64:
 1  4
 2  5

A is firstly a Matrix, then it becomes a adjoint. So A becomes a Union{Matrix{Float64}, Adjoint(Float64)} because it can be either of types/..

function eigSVD(A)
    transposed = false
    if size(A, 1) < size(A, 2)
        A = A'
        transposed = true
    end
    ....
end

This is (one of) the correct way to do so. We define a helper function eigenSVDHelper. Even input types to this funciton are different, the return types are always the same.

function eigenSVDHelper(A, transposed::Bool) 
    other code here...
    if transposed
        return V, S, U
    else
        return U, S, V
    end
end
function eigenSVD(A)
    transposed = false
    if size(A, 1) < size(A, 2)
         return eigenSVDHelper(A', true)
    else
         return eigenSVDHelper(A, false)
    end
end