3

我需要构建一个大的 NxN 稀疏带矩阵 A，其中 N = 570*720 = 410400（图像像素数）。

数学上，A(m,n) = C1 * exp(-|mn|^2); m = 1:N, n = 1:N

基本上它是在每一行评估的高斯函数，行索引是平均值和一些任意但小的标准偏差。

对于 N = 100，它看起来像：高斯

不幸的是，由于不必要的计算，它在 N = 410400 时运行非常慢。

1）使用for循环

A = sparse(N,N);
for i=1:N
A(i,:) = normpdf(1:N, i, 30);
end

由于调用 normpdf N 次，这是浪费且缓慢的。

2）计算 normpdf 一次 1:2N 平均在 N 然后循环移动 A 中的行与基于索引的适当平均值。matlab 中的 circshift 不能以不同的移位大小逐列移位矩阵。再次需要调用 circshift N 次 --> 浪费。

3）也许使用 mvnpdf 但它需要输入向量并且使用网格网格生成这些会
消耗大量内存。

4) 将 bsxfun 与用户定义的高斯函数（具有固定 SD）一起使用，因为 bsxfun 不接受 >3 个参数，因此接受两个参数。

关于如何有效实现这一点的任何想法？

谢谢

4

2 回答 2

4

我可以建议另一种方法。

什么会反对使用2*N-by-1可以通过对索引进行简单转换来索引的向量？像这样：

% Oli's solution, and your request: NxN matrix
N   = 100;
s   = 30;    
pdf = @(x) 1/(sqrt(2*pi)*s)*exp(-0.5*(bsxfun(@minus,x(:),1:N)/s).^2);
A   = pdf(1:N);

% My solution: 2*N x 1 vector
B = exp(-0.5*((-N:N)/s).^2) / s/sqrt(2*pi);

诀窍是找到一个很好的通用索引规则。这是如何做到的：

% Indexing goes like this:
fromB = @(ii,jj) B(N+1 + bsxfun(@minus, jj(:), ii)).';

ii = 30; 
jj = 23;

from_A = A(ii,jj)
from_B = fromB(ii,jj)

ii = 1:2;    
jj = 4:6;      

from_A = A(ii, jj)
from_B = fromB(ii,jj)

结果：

from_A =
   0.012940955690785
from_B =
   0.012940955690785

from_A =
   0.013231751582567   0.013180394696194   0.013114657203398
   0.013268557543798   0.013231751582567   0.013180394696194
from_B =
   0.013231751582567   0.013180394696194   0.013114657203398
   0.013268557543798   0.013231751582567   0.013180394696194

我相信您可以弄清楚如何进行冒号索引和使用end关键字之类的操作:)

于 2013-06-10T20:47:52.073 回答
3

首先，你真的不需要稀疏矩阵，除非你有至少 50% 的零但你的矩阵是full。

考虑一个普通的pdf

在此处输入图像描述

您可以轻松实现它，包括调用bsxfun()：

N   = 100;
s   = 30;
m   = 1:N;
pdf = @(x) 1/(sqrt(2*pi)*s)*exp(-0.5*(bsxfun(@minus,x(:),m)/s).^2);
B   = pdf(1:N);

一个简单的例子，mean = 1将sigma = 30阐明：

pdf = @(x) 1/(sqrt(2*pi)*30)*exp(-0.5*((x-1)/30).^2);
pdf(1)
ans =
        0.0132980760133811

normpdf(1, 1, 30)
ans =
        0.0132980760133811
于 2013-06-10T19:55:21.420 回答