Higher order SVD in Julia like Numpy?

General Usage
1

Hello all

I was wondering if there is any way to do higher order SVD in Julia using LinearAlgebra.jl like you can using numpy.linalg.svd or should I just use PyCall? I tried with LinearAlgebra.svd but it didn’t work out and I couldn’t find in the docs how to do it -

julia> A = randn(3, 5, 2)
3×5×2 Array{Float64, 3}:
[:, :, 1] =
  1.69017   0.539516  -0.436762  -0.882562   2.28308
 -0.180177  0.122244  -0.760002  -0.805284   1.15463
  0.218088  0.524262   1.12136   -2.33447   -0.706637

[:, :, 2] =
  0.850423  -1.12059    0.30392    -1.32743    -0.394143
  0.296605   0.766331  -0.25429     0.0635182   0.0745694
 -0.568156  -1.4385     0.0120621   1.66298     1.0844

julia> svd(A)
ERROR: MethodError: no method matching svd(::Array{Float64, 3})
Closest candidates are:
  svd(::StridedVecOrMat{T}; full, alg) where T<:Union{Float16, ComplexF16} at ~/.julia/juliaup/julia-1.8.0-rc3+0.x64/share/julia/stdlib/v1.8/LinearAlgebra/src/svd.jl:181
  svd(::StridedMatrix{T}, ::StridedMatrix{T}) where T<:Union{Float32, Float64, ComplexF32, ComplexF64} at ~/.julia/juliaup/julia-1.8.0-rc3+0.x64/share/julia/stdlib/v1.8/LinearAlgebra/src/svd.jl:407
  svd(::StridedVecOrMat{T}; full, alg) where T at ~/.julia/juliaup/julia-1.8.0-rc3+0.x64/share/julia/stdlib/v1.8/LinearAlgebra/src/svd.jl:178
  ...
Stacktrace:
 [1] top-level scope
   @ REPL[4]:1

Same thing in python -

In [1]: import numpy as np

In [2]: A = np.random.randn(3, 5, 2)

In [3]: np.linalg.svd(A)
Out[3]: 
(array([[[-0.6037455 ,  0.58416878,  0.2178343 , -0.36421969,
           0.33783197],
         [ 0.32082226, -0.42888505,  0.1985067 , -0.49056635,
           0.65808088],
         [ 0.25303468,  0.13119474,  0.92420027,  0.16216266,
          -0.19575117],
         [-0.4705167 , -0.38054641,  0.15731977,  0.655918  ,
           0.42287178],
         [ 0.49714378,  0.55926453, -0.18506123,  0.41250882,
           0.48544831]],
 
        [[-0.53278136,  0.63028902, -0.24075323,  0.4968516 ,
          -0.11855861],
         [-0.5110045 , -0.19041643, -0.59981223, -0.58222151,
           0.06212399],
         [-0.67157959, -0.32594425,  0.66462678,  0.00804075,
          -0.0307863 ],
         [-0.03712789, -0.67049536, -0.37229753,  0.63951442,
           0.03838423],
         [-0.05119373,  0.10331337,  0.04392   ,  0.07150708,
           0.98977943]],
 
        [[-0.50891794,  0.66685684,  0.31052877, -0.10016269,
           0.43571069],
         [-0.41478297, -0.66256106, -0.02431718, -0.43413115,
           0.44710934],
         [ 0.11320181, -0.25154376,  0.9476955 ,  0.07824908,
          -0.14021911],
         [-0.44217126, -0.22985906, -0.06724939,  0.86055279,
           0.08109132],
         [ 0.60052273, -0.01433001,  0.018202  ,  0.2341443 ,
           0.76420693]]]),
 array([[2.27220514, 0.61766628],
        [3.57499854, 1.38367195],
        [2.60157849, 2.01610413]]),
 array([[[ 0.42173425, -0.90671948],
         [-0.90671948, -0.42173425]],
 
        [[ 0.60493221, -0.79627697],
         [ 0.79627697,  0.60493221]],
 
        [[ 0.63651737, -0.77126237],
         [ 0.77126237,  0.63651737]]]))
1 Like
2

To be clear, this is just iterating the normal SVD for matrices over each slice of the array, right, not a generalization of SVD for higher-order tensors?

2 Likes
3

I need a generalization of SVD to higher order tensors. I apologize if I’m missing something, I am just starting to learn about higher order SVD and tensor networks, so I was unsure what numpy.linalg.svd was doing. I found ITensors.jl which seems to be doing what I require. It seems I got confused and just now realised numpy.linalg.svd was doing the iterative approach rather than a generalized SVD of higher-order tensors.

1 Like