Julia 或 R:通过最少的非对角线元素创建对称矩阵
Julia or R: create symmetric matrix by minimum of off-diagonal elements
这是一个复杂的问题。我在 R 中有一个任意方阵(也可以在 Julia 中),例如
> set.seed(420)
> A <- matrix(runif(16),nrow = 4,byrow = T)
> A
[,1] [,2] [,3] [,4]
[1,] 0.6055390 0.9702770 0.1744545 0.4757888
[2,] 0.7244812 0.8761027 0.3775037 0.6409362
[3,] 0.6546772 0.5062158 0.3033477 0.7162497
[4,] 0.2905202 0.1962252 0.3225786 0.8404279
我想把这个矩阵变换成对称矩阵,使得非对角线元素总是A
矩阵的2个对应的非对角线元素中的最小值。在上面的例子中,结果必须是:
[,1] [,2] [,3] [,4]
[1,] 0.6055390 0.7244812 0.1744545 0.2905202
[2,] 0.7244812 0.8761027 0.3775037 0.1962252
[3,] 0.1744545 0.3775037 0.3033477 0.3225786
[4,] 0.2905202 0.1962252 0.3225786 0.8404279
想要一种在 Julia 中高效且不耗时的编程方式 and/or R。必须适用于任何类型的方阵。
可能最好的方法是使用 Rcpp
进行几个 for 循环,但您也可以尝试一下:
idx <- A <= t(A)
A * idx + t(A) * (1L - idx)
编辑
基本的 R 解决方案在时间和内存方面都是低效的。如果您需要速度 and/or 内存,这里是使矩阵对称的 Rcpp
函数。一个就地修改它,另一个 returns 一个副本。
library(Rcpp)
cppFunction('
void inplaceSymmetric(NumericMatrix A) {
for (int i = 1; i < A.nrow(); ++i)
for (int j = 0; j < i; ++j)
if (A(i, j) < A(j, i))
A(j, i) = A(i, j);
else
A(i, j) = A(j, i);
}')
cppFunction('
NumericMatrix copySymmetric(NumericMatrix A) {
NumericMatrix C = clone(A);
for (int i = 1; i < A.nrow(); ++i)
for (int j = 0; j < i; ++j)
if (A(i, j) < A(j, i))
C(j, i) = C(i, j);
else
C(i, j) = C(j, i);
return C;
}')
在 Julia 中,一行是
B = [min(A[i,j], A[j,i]) for i in axes(A, 1), j in axes(A, 2)]
但是,这样做的工作量是所需的两倍。以下方法更有效,可以就地修改原始矩阵,如果您的矩阵不是正方形,则会生成一条很好的错误消息:
function symmetrize_min!(A::AbstractMatrix)
ax1 = axes(A, 1)
axes(A, 2) == ax1 || error("A must be square")
for j in ax1
for i in j+1:last(ax1)
A[i,j] = A[j,i] = min(A[i,j], A[j,i])
end
end
return A
end
在 Julia 中,修改参数以 !
结尾的函数是一种约定,作为对用户的警告。如果你不希望它修改 A
,那么你可以在调用它之前 copy(A)
或者做一个小的改变:
function symmetrize_min(A::AbstractMatrix)
ax1 = axes(A, 1)
axes(A, 2) == ax1 || error("A must be square")
B = similar(A)
for j in ax1
B[j,j] = A[j,j]
for i in j+1:last(ax1)
B[i,j] = B[j,i] = min(A[i,j], A[j,i])
end
end
return B
end
在 Julia 中,你应该拥抱循环——它们很快。而且因为它们很容易编写,所以看似复杂的问题变得简单。
这是一个复杂的问题。我在 R 中有一个任意方阵(也可以在 Julia 中),例如
> set.seed(420)
> A <- matrix(runif(16),nrow = 4,byrow = T)
> A
[,1] [,2] [,3] [,4]
[1,] 0.6055390 0.9702770 0.1744545 0.4757888
[2,] 0.7244812 0.8761027 0.3775037 0.6409362
[3,] 0.6546772 0.5062158 0.3033477 0.7162497
[4,] 0.2905202 0.1962252 0.3225786 0.8404279
我想把这个矩阵变换成对称矩阵,使得非对角线元素总是A
矩阵的2个对应的非对角线元素中的最小值。在上面的例子中,结果必须是:
[,1] [,2] [,3] [,4]
[1,] 0.6055390 0.7244812 0.1744545 0.2905202
[2,] 0.7244812 0.8761027 0.3775037 0.1962252
[3,] 0.1744545 0.3775037 0.3033477 0.3225786
[4,] 0.2905202 0.1962252 0.3225786 0.8404279
想要一种在 Julia 中高效且不耗时的编程方式 and/or R。必须适用于任何类型的方阵。
可能最好的方法是使用 Rcpp
进行几个 for 循环,但您也可以尝试一下:
idx <- A <= t(A)
A * idx + t(A) * (1L - idx)
编辑
基本的 R 解决方案在时间和内存方面都是低效的。如果您需要速度 and/or 内存,这里是使矩阵对称的 Rcpp
函数。一个就地修改它,另一个 returns 一个副本。
library(Rcpp)
cppFunction('
void inplaceSymmetric(NumericMatrix A) {
for (int i = 1; i < A.nrow(); ++i)
for (int j = 0; j < i; ++j)
if (A(i, j) < A(j, i))
A(j, i) = A(i, j);
else
A(i, j) = A(j, i);
}')
cppFunction('
NumericMatrix copySymmetric(NumericMatrix A) {
NumericMatrix C = clone(A);
for (int i = 1; i < A.nrow(); ++i)
for (int j = 0; j < i; ++j)
if (A(i, j) < A(j, i))
C(j, i) = C(i, j);
else
C(i, j) = C(j, i);
return C;
}')
在 Julia 中,一行是
B = [min(A[i,j], A[j,i]) for i in axes(A, 1), j in axes(A, 2)]
但是,这样做的工作量是所需的两倍。以下方法更有效,可以就地修改原始矩阵,如果您的矩阵不是正方形,则会生成一条很好的错误消息:
function symmetrize_min!(A::AbstractMatrix)
ax1 = axes(A, 1)
axes(A, 2) == ax1 || error("A must be square")
for j in ax1
for i in j+1:last(ax1)
A[i,j] = A[j,i] = min(A[i,j], A[j,i])
end
end
return A
end
在 Julia 中,修改参数以 !
结尾的函数是一种约定,作为对用户的警告。如果你不希望它修改 A
,那么你可以在调用它之前 copy(A)
或者做一个小的改变:
function symmetrize_min(A::AbstractMatrix)
ax1 = axes(A, 1)
axes(A, 2) == ax1 || error("A must be square")
B = similar(A)
for j in ax1
B[j,j] = A[j,j]
for i in j+1:last(ax1)
B[i,j] = B[j,i] = min(A[i,j], A[j,i])
end
end
return B
end
在 Julia 中,你应该拥抱循环——它们很快。而且因为它们很容易编写,所以看似复杂的问题变得简单。