forked from gherardovarando/gmat
-
Notifications
You must be signed in to change notification settings - Fork 1
/
qr.R
57 lines (40 loc) · 1.84 KB
/
qr.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
### implemented with QR decomposition
.rgmn_sqrt_qr <- function(rmsqrt, ug) {
p <- length(V(ug))
msqrt <- matrix(nrow = p, ncol = p, data = rmsqrt(p^2))
msqrt <- .ort_selective_qr(ug, msqrt)
return(msqrt %*% t(msqrt))
}
.ort_selective_qr <- function(ug, mort) {
madj <- igraph::as_adjacency_matrix(ug, type = "both", sparse = FALSE)
p <- nrow(madj)
degrees <- degree(ug)
order <- order(degrees, decreasing = FALSE)
# order <- 1:p
n_zeros <- length(degrees[degrees == 0])
### first we orthogonalize all the rows of the disconnected nodes (+1)
if (n_zeros > 0 ){
qrdec <- qr(t(mort[order[1:(n_zeros+1)],]))
mort[order[1:(n_zeros+1)],] <- t(qr.Q(qrdec))
}
## now we orthogonalize the other rows following the order
for (i in (n_zeros+2):p) {
row_iort <- which(madj[order[i], ] == FALSE)
row_iort <- row_iort[ row_iort %in% order[1:(i-1)]]
row_nort <- length(row_iort)
if (row_nort > 1) {
# mort[order[i], ] <- mort[order[i], ] -
# t( mort[row_iort,] ) %*% (solve( mort[row_iort,]%*% t(
# mort[row_iort,]) ) %*% mort[row_iort,] %*% mort[order[i],])
qrdec <- qr(t(mort[c(row_iort,order[i]),]))
mort[order[i],] <- t(qr.Q(qrdec))[row_nort + 1,]
# dec <- qr(t(mort[row_iort,]))
# Rr <- qr.R(dec)
# mort[order[i], ] <- mort[order[i], ] -
# t( mort[row_iort,] ) %*% (solve( t(Rr)%*% (Rr )) %*% mort[row_iort,] %*% mort[order[i],])
} else if (row_nort == 1) {
mort[order[i],] <- mort[order[i],] - mort[row_iort,]* sum(mort[order[i],]* mort[row_iort,])/ sum(mort[row_iort,]^2)
}
}
return (mort)
}