外部不在R工作

Question

我有3个向量x，y和z。我想创建一个矩阵，其第（i，j）个条目将是sum（（x [i] * z + y [j]）^ 2）

我正在尝试使用outer会发生错误。这里有一个例子：

x=y=1:10/10 
z=1:10 
outer(x,y,FUN=function(x,y,z) sum((x*z+y)^2),z) 

Error in outer(x, y, FUN = function(x, y, z) sum((x * z + y)^2), z) : 
    dims [product 100] do not match the length of object [1] 

Answer 1

您有以下总结：

\sum_{k=1}^n [(x_i*z_k + y_j)^2] 
\sum_{k=1}^n [(x_i*z_k)^2 + 2*x_i*z_k*y_j + y_j^2] 
\sum_{k=1}^n [(x_i*z_k)^2] + 2*x_i*y_j*\sum_{k=1}^n [z_k] + n*y_j^2 

因此，一个有效的方式来进行。将计算与​​x和z和计算之间的外积的第一项第二项与x和y之间的外部产品：

outer.sol <- function(x, y, z) { 
    n <- length(x) 
    xmult <- rowSums(outer(x, z, function(x, z) (x*z)^2)) 
    matrix(rep(xmult, n), nrow=n) + 
    outer(x, y, function(x, y) 2*sum(z)*x*y) + 
    matrix(rep(n*y^2, each=n), nrow=n) 
} 
outer.sol(x, y, z) 
#   [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] 
# [1,] 5.05 6.45 8.05 9.85 11.85 14.05 16.45 19.05 21.85 24.85 
# [2,] 17.70 20.20 22.90 25.80 28.90 32.20 35.70 39.40 43.30 47.40 
# [3,] 38.05 41.65 45.45 49.45 53.65 58.05 62.65 67.45 72.45 77.65 
# [4,] 66.10 70.80 75.70 80.80 86.10 91.60 97.30 103.20 109.30 115.60 
# [5,] 101.85 107.65 113.65 119.85 126.25 132.85 139.65 146.65 153.85 161.25 
# [6,] 145.30 152.20 159.30 166.60 174.10 181.80 189.70 197.80 206.10 214.60 
# [7,] 196.45 204.45 212.65 221.05 229.65 238.45 247.45 256.65 266.05 275.65 
# [8,] 255.30 264.40 273.70 283.20 292.90 302.80 312.90 323.20 333.70 344.40 
# [9,] 321.85 332.05 342.45 353.05 363.85 374.85 386.05 397.45 409.05 420.85 
# [10,] 396.10 407.40 418.90 430.60 442.50 454.60 466.90 479.40 492.10 505.00 

虽然你可以使用像mapply一个循环的解决方案，分别计算每一个元素，这将使你的代码一个很好的协议更容易阅读，注意，这会比完全量化的解决方案与outer慢得多：

mapply.sol <- function(x, y, z) { 
    n <- length(x) 
    matrix(mapply(function(i,j) sum((i*z+j)^2), rep(x, n), rep(y, each=n)), nrow=n) 
} 

# 1000 x 1000 case 
x <- y <- z <- 1:1000 
all.equal(outer.sol(x, y, z), mapply.sol(x, y, z)) 
# [1] TRUE 
system.time(outer.sol(x, y, z)) 
# user system elapsed 
# 0.060 0.001 0.060 
system.time(mapply.sol(x, y, z)) 
# user system elapsed 
# 18.522 1.986 20.432 

从评论它似乎真正的函数比问题中包含的二次函数更复杂，所以提供mapply.sol与您的函数是（相当缓慢）的方式。