-
Notifications
You must be signed in to change notification settings - Fork 15
/
svdDemo.R
87 lines (75 loc) · 2.4 KB
/
svdDemo.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
#' Demonstrate the SVD for a 3 x 3 matrix
#'
#' This function draws an \code{rgl} scene consisting of a representation of the identity matrix and a
#' 3 x 3 matrix \code{A}, together with the corresponding representation of the
#' matrices U, D, and V in the SVD decomposition,
#' A = U D V'.
#'
#' @param A A 3 x 3 numeric matrix
#' @param shape Basic shape used to represent the identity matrix: \code{"cube"} or \code{"sphere"}
#' @param alpha transparency value used to draw the shape
#' @param col Vector of 6 colors for the faces of the basic cube
#'
#' @return Nothing
#' @export
#' @author Original idea from Duncan Murdoch
#' @importFrom grDevices rainbow
#' @export
#' @examples
#' A <- matrix(c(1,2,0.1, 0.1,1,0.1, 0.1,0.1,0.5), 3,3)
#' svdDemo(A)
#'
#' \dontrun{
#' B <- matrix(c( 1, 0, 1, 0, 2, 0, 1, 0, 2), 3, 3)
#' svdDemo(B)
#'
#' # a positive, semi-definite matrix with eigenvalues 12, 6, 0
#' C <- matrix(c(7, 4, 1, 4, 4, 4, 1, 4, 7), 3, 3)
#' svdDemo(C)
#' }
#'
svdDemo <- function(A, shape=c("cube", "sphere"), alpha=0.7, col=rainbow(6)) {
d <- dim(A)
if( (length(d) != 2 ) | !all(d==3) ) stop("argument must be a 3 x 3 matrix")
svd <- svd(A)
U <- svd$u
D <- diag(svd$d)
V <- svd$v
shape <- match.arg(shape)
if (shape=="cube")
Shape <- cube3d(color=rep(col, each = 4), alpha=alpha)
else Shape <- subdivision3d(cube3d(color=rep(col, rep(4*4^4,6)), alpha=alpha), depth=4)
axes <- rbind( diag(3), -diag(3) )
rownames(axes) <- c("x", "y", "z", rep(" ", 3))
mult <- function(matrix, obj) transform3d(obj, t(matrix))
open3d()
# define the layout of subscenes
mat <- matrix(1:6, 2, 3)
layout3d(rbind(mat, mat+6), heights = c(3,1,3,1), sharedMouse=TRUE)
shade3d(Shape)
vectors3d(axes, frac.lab=1.2, headlength = 0.2, radius=1/20, lwd=3)
next3d()
text3d(0,0,0, "Identity", cex=1.5)
next3d()
next3d(reuse=FALSE)
next3d(reuse=FALSE)
shade3d(mult(A, Shape))
vectors3d(axes, frac.lab=1.2, headlength = 0.2, radius=1/20, lwd=3)
next3d()
text3d(0,0,0,"A", cex=1.5)
next3d()
shade3d(mult(U, Shape))
vectors3d(axes, frac.lab=1.2, headlength = 0.2, radius=1/20, lwd=3)
next3d()
text3d(0,0,0, "U", cex=1.5)
next3d()
shade3d(mult(D, Shape))
vectors3d(axes, frac.lab=1.2, headlength = 0.2, radius=1/20, lwd=3)
next3d()
text3d(0,0,0, "D", cex=1.5)
next3d()
shade3d(mult(V, Shape))
vectors3d(axes, frac.lab=1.2, headlength = 0.2, radius=1/20, lwd=3)
next3d()
text3d(0,0,0, "V", cex=1.5)
}