Lib for partial SVD computation

imorozov · September 21, 2020, 12:44pm

Hi,

Can someone recommend a FORTRAN library for partial SVD computation?
I.e. the case when only n first largest singular values are needed (corresponding U and VT are also needed).

I’ve tried DGESVD from lapack, but it seems there is no option to compute only several largest singular values or I’m missing something?

My DGESVD wrapper:

MODULE SVD
 USE, INTRINSIC :: ISO_C_BINDING,   ONLY: IK => C_INT, RK => C_DOUBLE, C_SIZEOF
  IMPLICIT NONE
  PRIVATE
  PUBLIC :: SVD_
  CONTAINS
  SUBROUTINE SVD_(NR, NC, MATRIX, SVD_LIST, U_MATRIX, V_MATRIX)
    INTEGER(IK), INTENT(IN) :: NR
    INTEGER(IK), INTENT(IN) :: NC
    REAL(RK), DIMENSION(NR, NC), INTENT(IN) :: MATRIX
    REAL(RK), DIMENSION(MIN(NR, NC)), INTENT(OUT) :: SVD_LIST
    REAL(RK), DIMENSION(NR, NR), INTENT(OUT) :: U_MATRIX
    REAL(RK), DIMENSION(NC, NC), INTENT(OUT) :: V_MATRIX
    REAL(8), DIMENSION(2*MAX(1, 3*MIN(INT(NR), INT(NC))+MAX(INT(NR), INT(NC)), 5*MIN(INT(NR), INT(NC)))) :: WORK
    INTEGER :: WORK_SIZE
    INTEGER :: INFO
    REAL(8), DIMENSION(NR, NC) :: COPY
    REAL(8), DIMENSION(MIN(NR, NC)) :: D_SVD_LIST
    REAL(8), DIMENSION(NR, NR) :: D_U_MATRIX
    REAL(8), DIMENSION(NC, NC) :: D_V_MATRIX
    COPY = REAL(MATRIX, 8)
    WORK_SIZE = SIZE(WORK)
    CALL DGESVD('A','A',INT(NR),INT(NC),COPY,INT(NR),D_SVD_LIST,D_U_MATRIX,INT(NR),D_V_MATRIX,INT(NC),WORK,WORK_SIZE,INFO)
    SVD_LIST = REAL(D_SVD_LIST, RK)
    U_MATRIX = REAL(D_U_MATRIX, RK)
    V_MATRIX = REAL(D_V_MATRIX, RK)
    V_MATRIX = TRANSPOSE(V_MATRIX)
  END SUBROUTINE SVD_
END MODULE SVD
PROGRAM EXAMPLE
  USE :: SVD
  IMPLICIT NONE
  INTEGER(IK) :: I
  ! SVD
  BLOCK
    REAL(RK) :: M(4_IK, 5_IK), S(4_IK), U(4_IK, 4_IK), V(5_IK, 5_IK)
    M = REAL(RESHAPE([1,4,7,8,2,5,8,9,3,6,9,10,4,7,10,11,5,8,11,12], SHAPE(M)), RK)
    WRITE(*, *) "INPUT MATRIX"
    DO I = 1_IK, 4_IK, 1_IK
      WRITE(*, *) M(I, :)
    END DO
    WRITE(*, *)
    CALL SVD_(4_IK, 5_IK, M, S, U, V)
    WRITE(*, *) "SINGULAR VALUES (SVD_)"
    WRITE(*, *) S
    WRITE(*, *)
    M = 0.0_RK
    DO I = 1_IK, SIZE(S, KIND = IK), 1_IK
      M(I, I) = S(I)
    END DO
    M = MATMUL(U, MATMUL(M, TRANSPOSE(V)))
    WRITE(*, *) "OUTPUT MATRIX"
    DO I = 1_IK, SIZE(S, KIND = IK), 1_IK
      WRITE(*, *) M(I, :)
    END DO
  END BLOCK
END PROGRAM EXAMPLE

lkedward · September 21, 2020, 1:07pm

I believe for truncated SVD you require an iterative procedure whereas LAPACK contains only direct algorithms (I think?).

You may want to look at ARPACK which provides iterative methods for eigenvalue and SVD problems. There’s an example here that shows “…how to use ARPACK to find a few of the largest singular values(sigma) and corresponding right singular vectors (v) for the the matrix A…”.

ARPACK uses a reverse communication interface to allow a matrix-free approach for very large sparse systems, so I would say it probably has a higher learning curve to get started than LAPACK though I haven’t yet had to use ARPACK myself.

ivanpribec · September 21, 2020, 1:24pm

There is SVDPACK: http://www.netlib.org/svdpack/
(I have no knowledge or personal experience with the interface)

The least squares package by Lawson & Hanson contains a SVD routine: https://www.netlib.org/lawson-hanson/all
Description of the algorithms are available in the book Solving Least Squares Problems.

This page describes two more options: SVDLIBC and PROPACK.

certik · September 21, 2020, 3:00pm

I use Arpack, here is a modern Fortran interface to the serial and parallel eigensolver:

I don’t know if it can do SVD.

imorozov · September 22, 2020, 2:09am

Thank you for your suggestions.
I’ve successfully tried to building ARPACK, it’s usage is quite cryptic.
Not much luck with other libs suggested by Ivan.

This explains it, indeed computation is full in LAPACK, good to know it.

ivanpribec · September 22, 2020, 9:54am

A few more libraries are mentioned here:

For the SVDLIBC you need to write a C interface. Is your matrix sparse?

imorozov · September 22, 2020, 11:20am

No, but in one of applications it is a Toeplitz matrix (transpose of it).
I’ve managed to rewrite ARPACK example and it should work for A(m,n) with m>n and requested number of singular values s<n, some polishing is still required, replace matmul with dgemv.

! gfortran -o dsvd -std=f2018 -Wall -pedantic -O3 -ffast-math -march=native -frecursive dsvd.f90 -lblas -larpack
! original code, https://github.com/opencollab/arpack-ng/blob/master/EXAMPLES/SVD/dsvd.f
! p110, http://li.mit.edu/Archive/Activities/Archive/CourseWork/Ju_Li/MITCourses/18.335/Doc/ARPACK/Lehoucq97.pdf
! to do, replace matmul() with dgemv (?), remove() transpose()
module svd

 use, intrinsic :: iso_c_binding, only: ik => c_int, rk => c_double
 implicit none

 private

 integer, public, parameter :: maxm = 512 ! -- max # of rows a(maxm,maxn)
 integer, public, parameter :: maxn = 256 ! -- max # of cols a(maxm,maxn)
 integer, public, parameter :: maxnev = 32 ! -- max # of singular values to compute (also max # of right singular vectors), maxnev <= maxn (?)
 integer, public, parameter :: maxncv = 16 ! -- largest # of basis vectors in in the Implicitly Restarted Arnoldi Process, optimal value (?)
 real(rk), public, parameter :: svd_level = 1.0E-6_rk ! -- singular value threshold
 integer, parameter :: ldu = maxm
 integer, parameter :: ldv = maxn

 public :: svd_

 external :: dsaupd
 external :: dseupd

 contains

 subroutine svd_(nr,nc,ns,matrix,list,rvec,lvec)
   integer, intent(in) :: nr ! number of rows
   integer, intent(in) :: nc ! number of cols
   integer, intent(in) :: ns ! number of singular values to compute
   real(rk), dimension(nr,nc), intent(in) :: matrix ! input matrix
   real(rk), dimension(ns), intent(out) :: list ! list of singular values
   real(rk), dimension(nc,ns), intent(out) :: rvec ! right singular vectors
   real(rk), dimension(nr,ns), intent(out) :: lvec ! left singular vectors
   real(rk), dimension(nc,nr) :: tmatrix ! transpose of input matrix
   real(rk), dimension(ns,ns) :: diag ! diagonal matrix
   integer :: i

   ! local arrays
   real(rk) :: v(ldv, maxncv) = 0.0_rk
   real(rk) :: workl(maxncv*(maxncv+8))
   real(rk) :: workd(3*maxn)
   real(rk) :: s(maxncv,2) = 0.0_rk
   real(rk) :: resid(maxn)
   real(rk) :: ax(maxm)
   logical :: select(maxncv)
   integer :: iparam(11)
   integer :: ipntr(11)

   ! local scalars
   character :: bmat*1
   character :: which*2
   integer :: ido
   integer :: m
   integer :: n
   integer :: nev
   integer :: ncv
   integer :: lworkl
   integer :: info
   integer :: ierr
   integer :: ishfts
   integer :: maxitr
   integer :: mode1
   integer :: nconv
   real(rk) :: tol
   real(rk) :: sigma

   ! transpose
   tmatrix = transpose(matrix)

   ! set dimensions
   m = nr
   n = nc

   ! arpack parameters
   nev   = ns ! # of singular values to compute, nev < n (?) nev == n (?)
   ncv   = n ! length of arnoldi factorization
   bmat  = 'I' ! OP = A'A
   which = 'LM' ! compute largest (magnitude) singular values, also LA

   ! arpack stopping rules and initials
   lworkl = ncv*(ncv+8) ! work array
   tol = 0.0_rk ! tolerance
   info = 0 ! initial error code
   ido = 0 ! reverse communication parameter

   ! arpack algorithm
   ishfts = 1   ! 0/1
   maxitr = 256 ! max number of iteraions
   mode1 = 1
   iparam(1) = ishfts
   iparam(3) = maxitr ! on out, actual number of iterations
   iparam(7) = mode1  ! mode

   ! main loop
   do
     call dsaupd(ido,bmat,n,which,nev,tol,resid,ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,info)
     if (.not.(ido .eq. -1 .or. ido .eq. 1)) exit
     call dot (m, n, matrix, workd(ipntr(1)), ax)
     call tdot (m, n, tmatrix, ax, workd(ipntr(2)))
   end do

   ! compute singular vectors
   call dseupd (.true.,'All',select,s,v,ldv,sigma,bmat,n,which,nev,tol,resid,ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,ierr)

   ! scale, note abs added for small negative vals
   nconv =  iparam(5)
   s(1:nconv,1) = sqrt(abs(s(1:nconv,1)))

   ! reverse
   s(1:nev,1) = s(nev:1:-1,1)

   ! set singular values
   list = s(1:nev,1)

   ! remove small
   where (list < svd_level) list = 0.0_rk

   ! set right vectors
   rvec(:,1:nev:1) = v(1:n,nev:1:-1)

   ! compute and set left vectors, u = a v s^-1
   diag = 0.0_rk
   do i = 1 , nev, 1
     if(list(i) > svd_level) diag(i, i) = 1.0_rk/list(i)
   end do
   lvec = matmul(matrix,matmul(rvec,diag))

 end subroutine svd_

 subroutine dot(row,col,matrix,x,w)
   integer, intent(in) :: row
   integer, intent(in) :: col
   real(rk), intent(in) :: matrix(row,col)
   real(rk), intent(in) :: x(col)
   real(rk), intent(out) :: w(row)
   w = matmul(matrix,x)
 end subroutine dot

 subroutine tdot(row,col,matrix,w,y)
   integer, intent(in) :: row
   integer, intent(in) :: col
   real(rk), intent(in) :: matrix(col,row)
   real(rk), intent(in) :: w(row)
   real(rk), intent(out) :: y(col)
   y = matmul(matrix,w)
 end subroutine tdot

end module svd

program test
  use, intrinsic :: iso_c_binding, only: ik => c_int, rk => c_double
  use svd
  implicit none

  integer, parameter :: nr = 5
  integer, parameter :: nc = 4
  integer, parameter :: ns = 3
  real(rk) ::  matrix(nr,nc)
  real(rk) ::  list(ns), diag(ns,ns)
  real(rk) ::  rvec(nc,ns)
  real(rk) ::  lvec(nr,ns)

  integer :: i

  matrix = 0.0_rk
  matrix(1,:) = real([1,2,3,4],rk)
  matrix(2,:) = real([2,3,4,5],rk)
  matrix(3,:) = real([3,4,5,6],rk)
  matrix(4,:) = real([4,5,6,7],rk)
  matrix(5,:) = real([5,6,7,8],rk)
  write(*, *) "input matrix"
  do i=1,5,1
    write(*, *) matrix(i,:)
  end do
  write(*, *)

  call svd_(nr,nc,ns,matrix,list,rvec,lvec)
  diag = 0.0_rk
  write(*, *) "singular values"
  do i = 1 , ns, 1
    diag(i, i) = list(i)
    write(*, *) diag(i,:)
  end do
  write(*,*)

  write(*,*) "right vectors"
  do i = 1, nc, 1
    write(*,*) rvec(i,:)
  end do
  write(*,*)

  write(*,*) "left vectors"
  do i = 1, nr, 1
    write(*,*) lvec(i,:)
  end do
  write(*,*)

  matrix = matmul(lvec,matmul(diag,transpose(rvec)))
  write(*,*) "U.S.V^T"
  do i=1, nr, 1
    write(*,*) matrix(i,:)
  end do
  write(*,*)
end program test

lkedward · September 22, 2020, 12:24pm

If you’re using gfortran you may save yourself some coding time by compiling with -fexternal-blas which will generate a call to BLAS for your matmul use. Then link against your preferred BLAS implementation as usual at link time.

imorozov · September 23, 2020, 6:38am

Great, interesting feature. But in this case, explicit calls should be better, since transpose of input matrix is used in multiplication.

Here is a version with BLAS matrix multiplication (tested for 2049x2048 matrices).
Left vector computed directly. Parameters SVD_BASIS and SVD_LENGH related to Arnoldi staff might not be optimal.

! gfortran -o dsvd -std=f2018 -Wall -pedantic -O3 -ffast-math -march=native -frecursive dsvd.f90 -lblas -larpack
! original code, https://github.com/opencollab/arpack-ng/blob/master/EXAMPLES/SVD/dsvd.f
! p110, http://li.mit.edu/Archive/Activities/Archive/CourseWork/Ju_Li/MITCourses/18.335/Doc/ARPACK/Lehoucq97.pdf

MODULE SVD
  USE, INTRINSIC :: ISO_C_BINDING, ONLY: IK => C_INT, RK => C_DOUBLE
  IMPLICIT NONE
  PRIVATE
  REAL(RK), PUBLIC, PARAMETER :: SVD_LEVEL = 1.0E-9_RK  ! SINGULAR VALUE THRESHOLD
  INTEGER,  PUBLIC, PARAMETER :: SVD_ROW = 2_IK**12_IK  ! MAX NUMBER OF ROWS
  INTEGER,  PUBLIC, PARAMETER :: SVD_COL = 2_IK**12_IK  ! MAX NUMBER OF ROWS
  INTEGER,  PUBLIC, PARAMETER :: SVD_BASIS = 64_IK      ! MAX NUMBER OF BASIS VECTORS IN THE IMPLICITLY RESTARTED ARNOLDI PROCESS, OPTIMAL VALUE (?)
  INTEGER,  PUBLIC, PARAMETER :: SVD_LENGH = 16_IK      ! LENGTH OF ARNOLDI FACTORIZATION
  INTEGER,  PUBLIC, PARAMETER :: SVD_LOOP = 256_IK      ! MAX NUMBER OF MAIN LOOP ITERAIONS
  PUBLIC :: SVD_TRUNCATED_
  PUBLIC :: MATRIX_
  EXTERNAL :: DSAUPD
  EXTERNAL :: DSEUPD
  EXTERNAL :: DGEMV
  EXTERNAL :: DGEMM
  CONTAINS
  ! TRUNCATED SVD (ARPACK)
  ! SVD_TRUNCATED_(<NR>,<NC>,<NS>,<MATRIX>(<NR>,<NC>),<LIST>(<NS>),<RVEC>(<NC>,<NR>),<LVEC>(<NR>,<NC>))
  SUBROUTINE SVD_TRUNCATED_(NR, NC, NS, MATRIX, LIST, RVEC, LVEC)
    INTEGER(IK), INTENT(IN) :: NR
    INTEGER(IK), INTENT(IN) :: NC
    INTEGER(IK), INTENT(IN) :: NS
    REAL(RK), DIMENSION(NR, NC), INTENT(IN) :: MATRIX
    REAL(RK), DIMENSION(NS), INTENT(OUT) :: LIST
    REAL(RK), DIMENSION(NC, NS), INTENT(OUT) :: RVEC
    REAL(RK), DIMENSION(NR, NS), INTENT(OUT) :: LVEC
    REAL(RK), DIMENSION(NS, NS) :: DIAG
    REAL(RK), DIMENSION(NC, NS) :: COPY
    INTEGER(IK) :: I
    REAL(RK), DIMENSION(SVD_COL, SVD_BASIS) :: V = 0.0_RK
    REAL(RK), DIMENSION(SVD_BASIS*(SVD_BASIS+8_IK)) :: WL
    REAL(RK), DIMENSION(3_IK*SVD_COL) :: WD
    REAL(RK), DIMENSION(SVD_BASIS,2_IK) :: S = 0.0_RK
    REAL(RK), DIMENSION(SVD_COL) :: ID
    REAL(RK), DIMENSION(SVD_ROW) :: AX
    INTEGER(IK) :: PAR(11_IK), PNT(11_IK)
    LOGICAL :: SELECT(SVD_BASIS)
    CHARACTER(LEN=1) :: MAT
    CHARACTER(LEN=2) :: MODE
    INTEGER(IK) :: IDO, NEV, NCV, WORK, INFO, IERR, NCONV
    REAL(RK) :: TOL, SIGMA
    ! ARPACK
    NEV    = NS                    ! # OF SINGULAR VALUES TO COMPUTE, NEV < N
    NCV    = MIN(SVD_LENGH, NC)    ! LENGTH OF ARNOLDI FACTORIZATION
    MAT    = 'I'                   ! OP = A^T.A
    MODE   = 'LM'                  ! COMPUTE LARGEST (MAGNITUDE) SINGULAR VALUES, ALSO LA
    WORK   = NCV*(NCV+8_IK)        ! WORK ARRAY SIZE
    TOL    = 0.0_RK                ! TOLERANCE
    INFO   = 0_IK                  ! INITIAL ERROR CODE
    IDO    = 0_IK                  ! REVERSE COMMUNICATION PARAMETER
    PAR(1) = 1_IK                  ! SHIFT, 0/1
    PAR(3) = SVD_LOOP              ! MAX NUMBER OF ITERAIONS
    PAR(7) = 1_IK                  ! MODE
    ! MAIN LOOP
    MAIN : DO
      CALL DSAUPD(IDO,MAT,NC,MODE,NEV,TOL,ID,NCV,V,SVD_COL,PAR,PNT,WD,WL,WORK,INFO)
      IF (.NOT.(IDO .EQ. -1_IK .OR. IDO .EQ. 1_IK)) EXIT
      CALL DGEMV('N',NR,NC,1.0_RK,MATRIX,NR,WD(PNT(1)),1_IK,0.0_RK,AX,1_IK)
      CALL DGEMV('T',NR,NC,1.0_RK,MATRIX,NR,AX,1_IK,0.0_RK,WD(PNT(2)),1_IK)
    END DO MAIN
    ! RIGHT SINGULAR VERCTORS
    CALL DSEUPD (.TRUE.,'ALL',SELECT,S,V,SVD_COL,SIGMA,MAT,NC,MODE,NEV,TOL,ID,NCV,V,SVD_COL,PAR,PNT,WD,WL,WORK,IERR)
    ! SCALE, REVERSE AND SET SINGULAR VALUES
    NCONV =  PAR(5)
    LIST = SQRT(ABS(S(NEV:1_IK:-1_IK,1_IK)))
    ! TRIM
    WHERE (LIST < SVD_LEVEL) LIST = 0.0_RK
    ! REVERSE AND SET RIGHT VECTORS
    RVEC(:,1_IK:NEV:1_IK) = V(1_IK:NC,NEV:1_IK:-1_IK)
    ! COMPUTE LEFT VECTOR, U = A.V.S^-1, DIRECT COMPUTATION
    DIAG = 0.0_RK
    DO I = 1_IK , NEV, 1_IK
      IF(LIST(I) > SVD_LEVEL) DIAG(I, I) = 1.0_RK/LIST(I)
    END DO
    CALL DGEMM('N','N',NC,NS,NS,1.0_RK,RVEC,NC,DIAG,NS,0.0_RK,COPY,NC)
    CALL DGEMM('N','N',NR,NS,NC,1.0_RK,MATRIX,NR,COPY,NC,0.0_RK,LVEC,NR)
  END SUBROUTINE SVD_TRUNCATED_
  SUBROUTINE MATRIX_(LENGTH, SEQUENCE, MATRIX)
    INTEGER(IK), INTENT(IN) :: LENGTH
    REAL(RK), DIMENSION(LENGTH), INTENT(IN) :: SEQUENCE
    REAL(RK), DIMENSION(LENGTH/2_IK+1_IK, LENGTH/2_IK), INTENT(OUT) :: MATRIX
    INTEGER(IK) :: I
    DO I = 1_IK, LENGTH/2_IK+1_IK, 1_IK
      MATRIX(I, :) = SEQUENCE(I:I-1_IK+LENGTH/2_IK)
    END DO
  END SUBROUTINE MATRIX_
END MODULE SVD

PROGRAM DSVD
  USE, INTRINSIC :: ISO_C_BINDING, ONLY: IK => C_INT, RK => C_DOUBLE
  USE SVD
  IMPLICIT NONE
  BLOCK
    INTEGER(IK), PARAMETER :: LENGTH = 2_IK**4
    INTEGER(IK), PARAMETER :: NR = LENGTH/2_IK + 1_IK
    INTEGER(IK), PARAMETER :: NC = LENGTH/2_IK
    INTEGER(IK), PARAMETER :: NS = 3_IK
    REAL(RK) :: SEQUENCE(LENGTH)
    REAL(RK) :: MATRIX(NR,NC)
    REAL(RK) :: SVD_LIST(NS)
    REAL(RK) :: DIAG(NS,NS)
    REAL(RK) :: RVEC(NC,NS)
    REAL(RK) :: LVEC(NR,NS)
    INTEGER :: I
    SEQUENCE = REAL([(I, I = 1, LENGTH, 1)], RK)
    CALL MATRIX_(LENGTH, SEQUENCE, MATRIX)
    WRITE(*, *) "INPUT MATRIX"
    WRITE(*, *) SIZE(MATRIX,1)
    WRITE(*, *) SIZE(MATRIX,2)
    DO I=1,NR,1
      WRITE(*, *) INT(MATRIX(I,:))
    END DO
    WRITE(*, *)
    CALL SVD_TRUNCATED_(NR, NC, NS, MATRIX, SVD_LIST, RVEC, LVEC)
    DIAG = 0.0_RK
    WRITE(*, *) "SINGULAR VALUES"
    DO I = 1, NS, 1
      DIAG(I, I) = SVD_LIST(I)
      WRITE(*, *) DIAG(I, :)
    END DO
    WRITE(*, *)
    WRITE(*, *) "RIGHT VECTORS"
    DO I = 1, NC, 1
      WRITE(*, *) RVEC(I,:)
    END DO
    WRITE(*, *)
    WRITE(*, *) "LEFT VECTORS"
    DO I = 1, NR, 1
      WRITE(*, *) LVEC(I, :)
    END DO
    WRITE(*, *)
    MATRIX = MATMUL(LVEC, MATMUL(DIAG, TRANSPOSE(RVEC)))
    WRITE(*, *) "U.S.V^T"
    DO I=1, NR, 1
      WRITE(*, *) CEILING(MATRIX(I, :))
    END DO
    WRITE(*, *)
  END BLOCK
END PROGRAM DSVD

lkedward · September 23, 2020, 7:55am

I believe -fexternal-blas can actually generate dgemm code for matmul with tranpose, see slide 15 in: Front-end optimization in gfortran, Thomas König, FortranCon 2020.

I’m glad to see you’ve got it working with ARPACK and thanks for sharing your code!

ivanpribec · September 23, 2020, 9:42am

I’ve found one more library called primme. A description of the Fortran interface can be found here.

Topic		Replies	Views
Singular Value Decomposition Implementation	6	1349	April 13, 2021
Optimal way to code getting spectra of matrices?	9	3091	September 1, 2021
Help with Blas band matrix operations (SOLVED) Help	10	1747	May 4, 2021
Issues with Pseudo-Inverse implementation using SVD	17	1279	June 21, 2023
LAPACK interfaces	45	3269	August 21, 2023

Lib for partial SVD computation

Related topics