<p>由于对这一点存在一些混淆,为了澄清问题,需要在注释末尾输入A的情况下找到p</p>
<p><strong>1)本征值</strong>假设A是对称正定矩阵,在问题中它是:</p>
<pre><code>P <- with(eigen(A), diag(1/sqrt(values)) %*% t(vectors))
# check
P %*% A %*% t(P)
## [1,] 1.000000e+00 -6.661338e-16 2.775558e-17
## [2,] -6.661338e-16 1.000000e+00 3.469447e-16
## [3,] 2.775558e-17 1.994932e-16 1.000000e+00
</code></pre>
<p><strong>2)optim</strong>另一种方法是数值求解这6个方程。请注意,参数化将P约束为对称,从而将参数数量减少到6。一个完整的3x3=9参数可能会导致数值问题,因为如注释中所述,如果P是一个解,那么HP对于任何正交H都太过重要</p>
<pre><code># input is upper triangle of P incl diagonal and output is P
p2P <- function(p) {
P <- diag(3)
P[upper.tri(P, diag = TRUE)] <- p
P + t(P) - diag(diag(P))
}
# sum of squares of difference between PAP' and diag(3)
f <- function(p) {
P <- p2P(p)
sum(c(P %*% A %*% t(P) - diag(3))^2)
}
res <- optim(1:6, f, method = "BFGS")
str(res)
## List of 5
## $ par : num [1:6] -0.2605 0.2573 0.5745 -0.0776 -0.027 ...
## $ value : num 4.69e-10
## $ counts : Named int [1:2] 167 49
## ..- attr(*, "names")= chr [1:2] "function" "gradient"
## $ convergence: int 0
## $ message : NULL
# check that PAP' is diag(3)
P <- p2P(res$par)
P %*% A %*% t(P)
## [,1] [,2] [,3]
## [1,] 9.999948e-01 -8.654260e-08 8.333054e-06
## [2,] -8.654260e-08 9.999955e-01 3.266114e-07
## [3,] 8.333054e-06 3.266114e-07 9.999832e-01
</code></pre>
<p><strong>3)nls</strong>我们可以交替使用<code>nls</code>来解方程。我们从上面使用p2P,但不需要f。请注意,P也是对称的</p>
<pre><code>result2 <- nls(c(diag(3)) ~ c(p2P(p) %*% A %*% t(p2P(p))),
start = list(p = 1:6), algorithm = "port")
result2
## Nonlinear regression model
## model: c(diag(3)) ~ c(p2P(p) %*% A %*% t(p2P(p)))
## data: parent.frame()
## p1 p2 p3 p4 p5 p6
## 0.26049 -0.25733 -0.57448 0.07757 0.02702 0.22665
## residual sum-of-squares: 1.354e-22
##
## Algorithm "port", convergence message: absolute function convergence (6)
# check that PAP' equals diag(3)
P <- p2P(coef(result2))
P %*% A %*% t(P)
## [,1] [,2] [,3]
## [1,] 1.000000e+00 7.954436e-13 5.084821e-14
## [2,] 7.954401e-13 1.000000e+00 -1.869643e-12
## [3,] 5.073719e-14 -1.869574e-12 1.000000e+00
</code></pre>
<p><strong>4)nls-2</strong>这也使用nls,但使用不同的p参数化,即对角矩阵(3个参数)乘以斜对称矩阵(3个参数)的指数。请注意,斜对称矩阵的指数是正交的,因此参数化将P约束为对角乘以正交</p>
<pre><code>library(expm)
p2P4 <- function(p) {
X <- matrix(0, 3, 3)
X[upper.tri(X)] <- p[1:3]
diag(p[4:6]) %*% expm(X - t(X))
}
res4 <- nls(c(diag(3)) ~ c(p2P4(p) %*% A %*% t(p2P4(p))), start = list(p = 1:6),
algorithm = "port")
res4
## Nonlinear regression model
## model: c(diag(3)) ~ c(p2P4(p) %*% A %*% t(p2P4(p)))
## data: parent.frame()
## p1 p2 p3 p4 p5 p6
## 3.3142 1.9743 1.6253 0.6500 0.1963 0.3663
## residual sum-of-squares: 2.863e-30
##
## Algorithm "port", convergence message: X-convergence (3)
# check that PAP' equals diag(3)
P <- p2P4(coef(res4))
P %*% A %*% t(P)
## [,1] [,2] [,3]
## [1,] 1.000000e+00 9.367507e-16 -2.827599e-16
## [2,] 9.159340e-16 1.000000e+00 -6.661338e-16
## [3,] -2.289835e-16 -6.106227e-16 1.000000e+00
</code></pre>
<h2>注</h2>
<p>我们将其用于:</p>
<pre><code>Lines <- "
10.166746 -2.619763 -6.717475
-2.619763 3.291888 3.052898
-6.717475 3.052898 22.313351"
A <- as.matrix(read.table(text = Lines))
</code></pre>