SVD projection, solve for U in X = UDV'

[neilcaithness]

I'm trying to project new data rows into a prior SVD space. ie. I want to solve for u in x = udv' which I think is u = xvd(-1) (e.g. from here https://blog.statsbot.co/singular-value-decomposition-tutorial-52c695315254).

From the output below I guess I must be wrong. Can some kind person help please?

x <- scale(iris[-5])

# svd
eg <- svd(x)
u <- eg$u
d <- diag(eg$d)
v <- eg$v

# check x = udv'
head(x)
#>      Sepal.Length Sepal.Width Petal.Length Petal.Width
#> [1,]   -0.8976739  1.01560199    -1.335752   -1.311052
#> [2,]   -1.1392005 -0.13153881    -1.335752   -1.311052
#> [3,]   -1.3807271  0.32731751    -1.392399   -1.311052
#> [4,]   -1.5014904  0.09788935    -1.279104   -1.311052
#> [5,]   -1.0184372  1.24503015    -1.335752   -1.311052
#> [6,]   -0.5353840  1.93331463    -1.165809   -1.048667
head(u %*% d %*% t(v))
#>            [,1]        [,2]      [,3]      [,4]
#> [1,] -0.8976739  1.01560199 -1.335752 -1.311052
#> [2,] -1.1392005 -0.13153881 -1.335752 -1.311052
#> [3,] -1.3807271  0.32731751 -1.392399 -1.311052
#> [4,] -1.5014904  0.09788935 -1.279104 -1.311052
#> [5,] -1.0184372  1.24503015 -1.335752 -1.311052
#> [6,] -0.5353840  1.93331463 -1.165809 -1.048667

# project onto u = xvd(-1)
head(x %*% v %*% d^-1)
#>      [,1] [,2] [,3] [,4]
#> [1,]  NaN  NaN  NaN  NaN
#> [2,]  Inf  NaN  NaN  NaN
#> [3,]  NaN  NaN  NaN  NaN
#> [4,]  NaN -Inf  NaN  NaN
#> [5,] -Inf -Inf -Inf -Inf
#> [6,]  NaN  NaN  NaN -Inf

^{Created on 2019-04-01 by the reprex package (v0.2.1)}

[Yarnabrina]

You're trying to find the inverse of the matrix d by d^(-1), which is wrong.

If you write A ^ k in R, where A is a matrix and k is a scalar, R will try to compute it elementwise, which obviously you don't want. Since d is diagonal, all non-diagonal entries will lead to Inf and hence you're facing this error.

To find the inverse, use the solve function. See below:

X <- scale(x = iris[-5])

svd_X <- svd(x = X)

U <- svd_X$u
D <- diag(x = svd_X$d)
V <- svd_X$v

all.equal(target = X,
          current = (U %*% D %*% t(x = V)),
          check.attributes = FALSE)
#> [1] TRUE

all.equal(target = U,
          current = (X %*% V %*% solve(a = D)))
#> [1] TRUE

^{Created on 2019-04-01 by the reprex package (v0.2.1)}

Hope this helps.

[neilcaithness]

Sorry, got it now. The matrix inverse is given by solve(d), not d^-1.

# project onto u = xvd(-1)
head(u)
#>             [,1]        [,2]         [,3]         [,4]
#> [1,] -0.10823953 -0.04099580  0.027218646  0.013710648
#> [2,] -0.09945776  0.05757315  0.050003401  0.058435855
#> [3,] -0.11299630  0.02920003 -0.009420891  0.016098333
#> [4,] -0.10989710  0.05101939 -0.019457133 -0.037416661
#> [5,] -0.11422046 -0.05524180 -0.003354363 -0.020379051
#> [6,] -0.09920300 -0.12718049 -0.005747892  0.003748828
head(x %*% v %*% solve(d))
#>             [,1]        [,2]         [,3]         [,4]
#> [1,] -0.10823953 -0.04099580  0.027218646  0.013710648
#> [2,] -0.09945776  0.05757315  0.050003401  0.058435855
#> [3,] -0.11299630  0.02920003 -0.009420891  0.016098333
#> [4,] -0.10989710  0.05101939 -0.019457133 -0.037416661
#> [5,] -0.11422046 -0.05524180 -0.003354363 -0.020379051
#> [6,] -0.09920300 -0.12718049 -0.005747892  0.003748828

^{Created on 2019-04-01 by the reprex package (v0.2.1)}

[neilcaithness]

Yes thanks very much. Just posted that when your answer came in.