显示标签为“科研笔记”的博文。显示所有博文
显示标签为“科研笔记”的博文。显示所有博文

2013-06-25

Machine Learning 第八波编程作业(完)——Anomaly Detection and Recommender Systems

仅列出核心代码:

1.estimateGuassian.m

mu = mean(X)';
X2 = (X - ones(m, 1)*mu').^2;
sigma2 = mean(X2);

2.selectThreshold.m

cvPredictions = (pval < epsilon);
tp = sum((cvPredictions == 1) & (yval == 1));
fp = sum((cvPredictions == 1) & (yval == 0));
fn = sum((cvPredictions == 0) & (yval == 1));
prec = tp/(tp + fp);
rec = tp/(tp + fn);
F1 = 2*prec*rec/(prec + rec);

3.cofiCostFunc.m

X1 = (X*Theta'- Y).*R;
reg1 = (sum(sum(X.^2)) + sum(sum(Theta.^2)))*lambda/2;
J = sum(sum((X1).^2))/2 + reg1;

X_grad = X1*Theta + lambda*X;
Theta_grad = X1'*X + lambda*Theta;

课程地址:https://www.coursera.org/course/ml

 

2013-06-18

Machine Learning 第七波编程作业——K-means Clustering and Principal Component Analysis

仅列出核心代码:

1.findClosestCentroids.m

m = size(X, 1);
len = zeros(K, 1);
for i = 1:m
    for j = 1:K
        len(j) = norm(X(i, :) - centroids(j, :))^2;
    end
    [~, idx(i)] = min(len);
end

2.computeCentroids.m

for k = 1:K
    ind = find(idx == k);
    centroids(k, :) = mean(X(ind, :));
end

3.pca.m

Sigma = X'*X/m;
[U,S,~] = svd(Sigma);

4.projectData.m

Z = X * U(:, 1:K);

5.recoverData.m

X_rec = Z * U(:, 1:K)';

课程地址:https://www.coursera.org/course/ml

2013-06-11

Machine Learning 第六波编程作业——Support Vector Machines

仅列出核心代码:

1.gaussianKernel.m

sim = exp(-sum((x1 - x2).^2) /(2*sigma^2));

2.dataset3Params.m

TD =  [0.01, 0.03, 0.1, 0.3, 1, 3, 10, 30];
pre_err = zeros(length(TD));
for i = 1:length(TD)
    for j = 1:length(TD)
        C = TD(i);
        sigma = TD(j);
        model= svmTrain(X, y, C, @(x1, x2) gaussianKernel(x1, x2, sigma));
        predictions = svmPredict(model, Xval);
        pre_err(i, j) = mean(double(predictions ~= yval));
    end
end
mm = min(min(pre_err));
[ind_C, ind_sigma] = find(pre_err == mm);
C = TD(ind_C);
sigma = TD(ind_sigma);

3.processEmail.m

for i = 1:length(vocabList)
v = strcmp(str, vocabList(i));
    if v==1
        word_indices = [word_indices ; i];
    end
end

4.emailFeatures.m

x(word_indices) = 1;

课程地址:https://www.coursera.org/course/ml

2013-06-03

Machine Learning 第五波编程作业 – Regularized Linear Regression and Bias/Variance

仅列出核心代码:

1.linearRegCostFunction.m

h = X * theta;
J = (X * theta - y).' * (X * theta - y) / (2*m)...
    +(lambda/(2*m)) * sum(theta(2:end).^2);
grad = grad(:);
grad(1) = (X(:, 1).' * (h - y)) /m;

grad(2:end) = (X(:, 2:end).' * (h - y)) /m ...
+ (lambda/m) * theta(2:end);

2.learningCurve.m

for i = 1:m
    Xi = X(1:i, :);
    yi = y(1:i);
    lambda = 1;
    [theta] = trainLinearReg(Xi, yi, lambda);
    lambda = 0;
    % For train error, make sure you compute it on the training subset
    [error_train(i), ~] = linearRegCostFunction(Xi, yi, theta, lambda);
    % For validation error, compute it over the entire cross validation set
    [error_val(i), ~] = linearRegCostFunction(Xval, yval, theta, lambda);
end

3.polyFeatures.m

for i =1:p
    X_poly(:, i) = X(:, 1).^i;
end

4.validationCurve.m

for i = 1:length(lambda_vec)
    [theta] = trainLinearReg(X, y, lambda_vec(i));
    % For train error, make sure you compute it on the training subset
    [error_train(i), ~] = linearRegCostFunction(X, y, theta, 0);
    % For validation error, compute it over the entire cross validation set
    [error_val(i), ~] = linearRegCostFunction(Xval, yval, theta, 0);
end

课程地址:https://www.coursera.org/course/ml

2013-05-28

Machine Learning 第四波编程作业 - Neural Networks: Learning

仅列出核心代码:

1.sigmoidGradient.m

h = 1.0 ./ (1.0 + exp(-z));
g = h.*(1 - h);

2.randInitializeWeights.m

epsilon_init = 0.12;
W = rand(L_out, 1 + L_in)*2*epsilon_init - epsilon_init;

3.nnCostFunction.m

% cost function

A1 = X;
A1 = [ones(m, 1), A1];
Z2 = A1 * Theta1.';
A2 = sigmoid(Z2);
A2 = [ones(m, 1), A2];
Z3 = A2 * Theta2.';
A3 = sigmoid(Z3);
H = A3;
Y = zeros(m, num_labels);
for ind = 1:m
Y(ind, y(ind)) = 1;
end
K = num_labels;
Jk = zeros(K, 1);

for k =1:K

Jk(k) = ( -Y(:, k).' *log(H(:, k)) )-( (1 - Y(:, k)).' * log(1-H(:, k)) );

end
J = sum(Jk)/m;

J = J + ( lambda/(2*m) )*( sum(sum(Theta1(:, 2:end).^2))+sum(sum(Theta2(:, 2:end).^2)) );

% Unroll gradients

delta3 = A3 - Y;
delta2 = delta3*Theta2.* (A2.*(1-A2));
delta2 = delta2(:, 2: end);
Delta2 = zeros(size(delta3, 2), size(A2, 2));
Delta1 = zeros(size(delta2, 2), size(A1, 2));
for i=1:m
Delta2 = Delta2 + delta3(i, :).' * A2(i, :);
Delta1 = Delta1 + delta2(i, :).' * A1(i, :);
end
Theta1_grad = Delta1/m;
Theta1_grad(:, 2:end) = Theta1_grad(:, 2:end) + Theta1(:, 2:end)*(lambda/m);
Theta2_grad = Delta2/m;
Theta2_grad(:, 2:end) = Theta2_grad(:, 2:end) + Theta2(:, 2:end)*(lambda/m);
grad = [Theta1_grad(:) ; Theta2_grad(:)];
end

课程地址:https://www.coursera.org/course/ml

2013-05-20

Machine Learning 第三波编程作业 – Multi-class Classification and Neural Networks

仅列出核心代码:

1.lrCostFunction.m

h = sigmoid(X * theta);   %   h_theta(X) : m*1
%   Cost func
J = (-log(h.')*y - log(ones(1, m) - h.')*(ones(m, 1) - y)) / m ...
    +(lambda/(2*m)) * sum(theta(2:end).^2);

%   Gradient
grad(1) = (X(:, 1).' * (h - y)) /m;

grad(2:end) = (X(:, 2:end).' * (h - y)) /m ...
    + (lambda/m) * theta(2:end);

2.oneVsAll.m

options = optimset('GradObj', 'on', 'MaxIter', 50);
initial_theta = zeros(size(X, 2), 1);
for c = 1:num_labels
    [all_theta(c, :)] = fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)),...
        initial_theta, options);   
end

3.predictOneVsAll.m

[~, p] = max(X * all_theta.', [], 2);

4.predict.m

X = [ones(size(X), 1), X]; % Add ones to the X data matrix

X1 = sigmoid(X * Theta1.');
X1 = [ones(size(X1), 1), X1]; % Add ones to the X1 data matrix

[~, p] = max(X1 * Theta2.', [], 2);

课程地址:https://www.coursera.org/course/ml

Machine Learning 第二波编程作业 – Logistic Regression

仅列出核心代码:

1.plotData.m

ind1 = find(y==1); ind0 = find(y==0);
plot(X(ind1, 1), X(ind1, 2), 'k+','LineWidth', 2, 'MarkerSize', 7);
plot(X(ind0, 1), X(ind0, 2), 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 7);

2.sigmoid.m

g = 1 ./ (ones(size(z)) + exp(-z));

3.costFunction.m

h = sigmoid(X * theta); % h_theta(X) : m*1
J = (-log(h.')*y - log(ones(1, m) - h.')*(ones(m, 1) - y)) / m;
grad = (X.' * (h - y)) /m;

4.predict.m

h = sigmoid(X * theta);
p = (h >= 0.5);

5.costFunctionReg.m

h = sigmoid(X * theta); % h_theta(X) : m*1
% Cost func
J = (-log(h.')*y - log(ones(1, m) - h.')*(ones(m, 1) - y)) / m ...
+(lambda/(2*m)) * sum(theta(2:end).^2);

% Gradient
grad(1) = (X(:, 1).' * (h - y)) /m;

grad(2:end) = (X(:, 2:end).' * (h - y)) /m ...
+ (lambda/m) * theta(2:end);


课程地址:https://www.coursera.org/course/ml

2013-05-10

Machine Learning 第一波编程作业 - Linear Regression

仅列出核心代码:

1.computeCost

J = sum((X*theta-y).^2)/(2*m);

2.gradientDescent

theta = theta - (1/m)*alpha*(X.'*(X*theta-y));

3.featureNormalize

mu = mean(X);
sigma = std(X);
X_norm = (X - ones(size(X, 1), 1) * mu) ./ (ones(size(X, 1), 1) * sigma);

4.computeCostMulti

J = (X * theta - y).' * (X * theta - y) / (2*m);

5.gradientDescentMulti

theta = theta - (1/m)*alpha*(X.'*(X*theta-y));

6.normalEqn

theta = inv(X.' * X) * X.' * y;

课程地址:https://www.coursera.org/course/ml

2013-05-07

《A Beginner's Guide to Irrational Behavior》课程论文

China is a labor-populous country, known as the "factory of the world". Some enterprises e.g. the Apple Intel etc., all set OEM factories located in China. I live in Chengdu, Sichuan in China. There is a huge factory with 110,000 workers, Foxconn. Two years ago, employee suicided in Foxconn Shenzhen factory frequently. During a short period(four months), twelve workers jumped from the top of building. As in such a tragic way to end their lives, it shocked the society. Soon, the news attracted attention of the media of China and abroad.

What's the reason of Foxconn's tragedy? Many experts (most of them were psychologists) had made a deep analysis. Some said that Foxconn put too much pressure onto its workers and the cost of living in such a big city in China is high. So that made a very obviously contrast between the two incentives. As a result of such effect, the workers' emotion inevitably lead to go extremely. On the other hand, there were some reports to disclose Foxconn industrial park's management, which was very mechanical, totally inhumanness. Among workers, there were also competitive relationships. It is difficult to establish friendship, unable to find a sense of belonging. Therefore, over time, those workers will be more and more unhappy....

In my opinion, one thing is certain: the workers were hard to feel a sense of accomplishment or happiness. Dan Ariely etc. (D Ariely, 2008) and Michael Norton etc. (M Norton, 2012), made a statement of a unanimous conclusion in their papers, which is: "a sense of identity" is important to everyone. This sense of identity comes from both the completion of the work and the affirmation of the results of the work from others. As "IKEA Effect" indicated, a doubly cherish would be put on the fruit of the task which need someone tried hard to achieve. And such endeavor would better be "systematically", for example, a furniture assembled by someone's own power. But if it is over-specialized, such as the division of labor(A Smith, 1776), there wouldn't be the similar effect. Foxconn's employees just like the products of socialization which makes the extreme division of labor. They are like a big car bed's bearing, working repeatedly  for times. They do not know what is the specific role of such assembly work, day after day, year after year. Where is the "meaning"? Whether it is to assemble the the Legos or assemble iPhone, people must be able to genuineness of seeing the fruits of their own labor. And that will take the initiative to stimulate a strong desire to finish their jobs. Otherwise, do not say to love the work, maybe they will feel tired of living for a long term.

Of course, as complicated "toys" as the iPhone, it is impossible to expect one person to assemble it correctly, and also the efficiency will be too low. But we can bring a sense of accomplishment to workers through other ways. For instance, in each specific part of the assembly, set a visualization object identifies which part of the mobile phone the workers' job belong to. Or in every step-assembly, set a graphical phone showing more completely, for every workers, a virtual mobile phone was assembled after several hours. And then, they would touched the goal of completion and feel some kind meaning. The principle of such mechanism is similar with the electric toothbrush's handle with smiley -- although it seems a little bit silly, sometimes it's very effective.

References:

Ariely, D., Kamenica, E. & Prelec, D. (2008). Man’s Search for Meaning: The Case of Legos. Journal of Economic Behavior and Organization, 67, 671-677.

Norton, M. I., Mochon, D., & Ariely, D. The IKEA Effect: When Labor Leads to Love. Harvard Business School Marketing Unit Working Paper, (11-091).

Smith A. The Wealth of Nations (1776)[J]. New York: Modern Library, 1937, 740.



最后得了7.5分。
一个比较长的评语的是:

I commend you on your use of English, which is generally clear in terms of meaning. You describe a problem of worker despair, leading to suicide. In paragraph 2, you refer to "a very obviously contrast between the two incentives." A clear statement of the two incentives and related behavior would be helpful. In addition, your assertion that inhumane treatment and lack of meaning in work were causative factors in the suicides. Your argument would be much stronger if you gave examples to work conditions to support these assertions.

想想我自己给别人打得分也都很高,瞬间淡定了。

2013-01-16

月球为何不是奶酪做的——再谈江本胜们何以如此流行

最近,在微博上看到有人用草莓重复了前一段风靡一时的「米饭行为实验」——就是对蒸熟的米饭说不同的词语(赞扬或者咒骂),数周之后,受到咒骂的米饭就会发霉变质,而得到赞扬的米饭则没有什么变化。这个实验的结果是令人震惊的——人的情感居然会影响物的质变。广州还有个什么「青少年危机干预中心」的张庆主任表示,「米饭实验得到了心理学和量子物理学的理论支持,是真实可信的。」(http://news.nfmedia.com/nfdsb/content/2011-11/30/content_34078140.htm)小朋友们做完这个实验也纷纷表示,以后要多多使用积极的词语,要养成一颗感恩的心。

这些看似有理有据,从实践得出的「真知」让我想起了一个日本人——江本胜。

江本胜是干嘛的?维基百科上称之为作家、商人,而众多书城则以「医学博士」的头衔来称呼他。12年前,他写了一本叫做《水知道答案》的小书,然后一炮而红。这本书由大量的照片组成,每一幅都是清晰的冰晶显微照片。江本胜称水在结晶的过程中会受到人类情感的影响,譬如:对其恶语相向,则结晶会变得丑陋;对其大加褒扬,则结晶分外好看;又或者,对其演奏优雅的古典音乐,则结晶亦显精致;如果对其播放噪声,则结晶一无是处……这些结论都是「有图有真相」,令人叹为观止。

不过且慢,这种看上去像模像样的「科学实验」靠谱么?

我想任何一个有着基本科学素养的人都会对这个匪夷所思的结论打一个大大的问号吧。

其实稍微想一想,就知道这个结论是站不住脚的。我从两方面来论证,先来复杂的:首先江本胜的实验不符合基本的科学实验要求,比如说没有设置对照组,没有引入双盲机制,更没有严格限定前提条件;其次江本胜的实验过程看似有理,实则幼稚,比如好话和坏话——这个是存在人的主观判断的,一个人的标准未必是另一个的标准,并且他仅仅实验了日语(或者常见的国际语言),那么那些生僻的语言呢,你是不是应该把地球上所有的语言都试验一遍才算严谨?最后江本胜的结论也给的非常牵强,冰晶的好看与否也是和人的主观审美情趣息息相关的,我就觉得几幅被噪声污染的冰晶挺好看,说不定也会有人觉得那些沐浴在赞扬声中的冰晶看上去显得丑陋。这些仅仅凭一厢情愿得出的结论,一无数据支撑、二无标准参照,简直漏洞百出。
(具体请参见:
http://www.guokr.com/article/60265/
还有这里:
http://discovery.163.com/special/waterisfool/
以及这里:
http://blog.sciencenet.cn/blog-1557-3412.html

那么接下来是简单的论证:Steve Mirsky在他的短文《Moon Not Made of Cheese, Physicist Explains》(http://www.scientificamerican.com/podcast/episode.cfm?id=moon-not-made-of-cheese-physicist-e-11-10-19)中引用了一个例子:「月亮为何不是由绿色奶酪做成的?」物理学家给出了一个简洁有力的回答:「因为这个想法是荒谬的。」得出这个结论的底气来自于「可以用我们对其他的事物的认知判断一个假设的合理性。」类似「月亮为何不是奶酪」、「人类为何没有任督二脉」、「占星术为何是一派胡言」、「谩骂为何会使得米饭变臭」等等问题,统统可以用「荒谬」二字来答复。

这是不是有些专断?人们喜欢说没有调查就没有发言权。那么我并没有亲自去蒸两碗米饭来检验一番,何以如此信誓旦旦地说它是荒谬的呢?

靠的是现有的科学理论和合理的逻辑推测。我无需亲自去重复实验,只需找出他实验中的漏洞,证明他的论证不靠谱,就完成了我的反驳。这是一个再简单不过的道理,可是总是有人拎不清。包括以科研素养著称的博士生群体,我曾经在课堂上与多位博士生辩论,试图告诉他们江本胜的书就是伪科学,可令我吃惊的是,有众多同学辩称「你又没有亲自验证怎么敢轻易下结论?」当我指出江本胜设计的实验根本站不住脚时,居然有人说「那是你不敢于跳出现有理论的框框。」我还能说什么呢……后来在与这些同学的线上交流中,我惊异地发现,他们有笃信什么「宇宙精神」的、有坚信特异功能的、有亲自练习硬气功的,令人大开眼界。面对着论文发了一篇又一篇可是基本的理性思维都不具备的高材生们,我无言以对。

所以,每次看到我们的媒体和出版者不负责任地把这些异想天开的劳什子包装一番然后印上「科普」的字眼放在书店的「科学」类书架内,我就无名火起三丈——这不是赤裸裸的欺诈么。不要说广大未经受严格科学训练的家长,即使他们经受了正规的科学训练又怎样?我们中的很多人已经失去了质疑的能力,失去了独立思考的习惯,并且一厢情愿地喜欢听那些「好听」的论点(譬如《水知道答案》中劝人向善的告诫),拿回去再一本正经地教给自己的小孩,然后小朋友们就会照猫画虎,去重复「米饭实验」……这样的「科普」,不是把人给带到沟里去了?

最后,我想说,尽管我们有了诸如果壳网、松鼠会这样活跃的科普团体,但我们依旧停留在向大众转发科学技巧的浅层科普,什么时候达到了卡尔·萨根那种把科学精神深深植入作品的地步,公众的科学素养才会有真正的提高。

米饭实验

2012-09-18

关于 Mathematica 的 Manipulate 函数无法释放CPU的解决办法

最近遇到一个小问题,就是在使用 Manipulate 函数时,一旦涉及到迭代、循环语句,代码就会反复执行(右侧方框加黑——而换用 Animate 函数则没有这种情况),只能手动中止,例如:
   
    Manipulate[
     date = Table[0, {num}];
     For[i = 1, i <= num, i++, date[[i]] = Sin[i]];
     ListPlot[{date}]
     , {{num, 2, "iteration num"}, 1, 20, 1}]

再来一个双层循环的例子:

    Manipulate[           
      buf = Table[0, {10}];
      Mw = Table[0, {10}];
      For[j = 1, j <= 10, j++,
       For[i = 1, i <= 10, i++,
        w = RandomVariate[NormalDistribution[0, var], 10];
        buf[[i]] = Mean[w];];
       Mw[[j]] = Mean@Cos[buf];];     
      ListLinePlot[Mw, PlotRange -> Full, DataRange -> All, Mesh -> Full,
       GridLines -> Automatic,
       GridLinesStyle -> Directive[Orange, Dashed], Frame -> True]     
     ,{{var, 3, "Var"}, 1, 5, 0.5}]
 
这两个例子还不是很严重,但是结构一旦再复杂一点, 就不太好办了,有时甚至会将前端卡死……

去网上搜索、发问,收到了一些具有启发性的回答。有人说到 Manipulate 是实时计算的,而 Animate 是统一计算所有值后再呈现结果,二者不同。

这样一来, 一旦 Manipulate 涉及到迭代、循环一类的语句,就会反复调用,无法退出,这也是逻辑上讲得通的,至少不使用迭代语句,Manipulate 就不会发生这种情况。

如果这是问题的原因,那么如何解决呢?

首先是修改算法,对于简单的问题,将循环语句转化为非循环语句或许行得通,可是对于很大一类迭代算法来说,这不是件轻易可以办到的事情,还是要从工具本身去需找答案。

于是我再去查看 Help文档,终于发现 Mathematica 的开发者早已考虑到这个问题,而且我遇到的问题也不能算作缺陷,只是一种非预设的情况罢了。

以下摘抄自「高级操作(Manipulate)功能」:

    每当步长被从零移开时,内容区会就会连续地更新,一个 CPU 监视器会指示 Mathematica 正在使用 CPU 时间. 你让这种情况持续多久,它就会持续多久.

    然而在某些情况下,持续的再计算是没有意义和不受欢迎的.
   
        例1:   
        Manipulate[
         temp = n;
         temp = temp^3;
         Graphics[{Thickness[0.01], Line[{{0, 0}, {n, temp}}]},
          PlotRange -> 1],
         {n, -1, 1}]
       
        例2:
        Manipulate[
         f[x_] := x^3;
         Graphics[{Thickness[0.01], Line[{{0, 0}, {n, f[n]}}]},
          PlotRange -> 1],
         {n, -1, 1}]
   
    在这两个情况下,这个问题都可以通过把引起问题的那些变量在一个 Module 里设成局部变量来解决. (这无论如何这都是一个好的编程习惯,远不止是为了避免无意义的更新)
   
        例1:
        Manipulate[Module[{temp},
          temp = n;
          temp = temp^3;
          Graphics[{Thickness[0.01], Line[{{0, 0}, {n, temp}}]},
           PlotRange -> 1]],
         {n, -1, 1}]
        
        例2:
        Manipulate[Module[{f},
          f[x_] := x^3;
          Graphics[{Thickness[0.01], Line[{{0, 0}, {n, f[n]}}]},
           PlotRange -> 1]],
         {n, -1, 1}]
    
    不管你对局部 Module 变量做什么都不会造成重新触发,因为这是 Module 定义的一部分,即一次调用后的变量值不会存留到下一次(所以下一轮的结果不会仅仅因为当前一轮运行中对局部变量所做的任何动作而有任何不同).
   
    另一个解决的办法是使用 Manipulate 的 TrackedSymbols (跟踪的符号)选项来控制哪一个变量能被允许引起更新行为. 默认的值, Full (全部),意味着所有在第一个自变量中明确出现(词汇的)的符号都会被跟踪. (这意味着,除了其它事情之外,在你使用的 Manipulate 例子中函数定义里的临时变量和其它如此的问题将不会引起无限重复计算问题,这是因为它们不在第一个自变量中明显地出现,而只是通过你调用的函数间接地出现.)

    来看第二个例子,如果由于某种原因你不想要 f 成为一个局部的 Module 变量,并且你不能把它的定义移到 Manipulate 之外(在更复杂的例子中,这两种情况有些时候都是很可能会出现的),你能用 TrackedSymbols 来取消由 f 触发的更新:

        Manipulate[
         f[x_] := x^3;
         Graphics[{Thickness[0.01], Line[{{0, 0}, {n, f[n]}}]},
          PlotRange -> 1],
         {n, -1, 1}, TrackedSymbols :> {n}]
    
    这个例子只在移动滑块从而改变 n 的数值时才更新内容区域.

   

以上内容完全解决了我的疑问,我将代码稍作改动,问题搞定,如下:

    Manipulate[Module[{i},
      date = Table[0, {num}];
      For[i = 1, i <= num, i++, date[[i]] = Sin[i]];
      ListPlot[{date}]]
     , {{num, 4, "iteration num"}, 1, 20, 1}]
    
    Manipulate[           
     buf = Table[0, {10}];
     Mw = Table[0, {10}];
     For[j = 1, j <= 10, j++,
      For[i = 1, i <= 10, i++,
       w = RandomVariate[NormalDistribution[0, var], 10];
       buf[[i]] = Mean[w];];
      Mw[[j]] = Mean@Cos[buf];];
     ListLinePlot[Mw, PlotRange -> Full, DataRange -> All, Mesh -> Full,
      GridLines -> Automatic, GridLinesStyle -> Directive[Orange, Dashed],
       Frame -> True]
     , {{var, 3, "Var"}, 1, 5, 0.5}
     , TrackedSymbols :> {var}]

正如高手所说的,工具本身的 Help 是最棒的专家系统。

此外,文档这一节的末尾提到这么一段话:「具体在什么时候一个给定的动态表达式会被更新这个话题是很复杂的,这在 "动态简介" 和 "高级动态功能" 中被提到了. 在阅读那些文件时,始终注意 Manipulate 只是在 Dynamic 中把它的第一个自变量包围起来并且把它的 TrackedSymbols 选项的数值传送个在那里面的 Refresh. 所有与更新有关的事情都是由那个 Dynamic 和 Refresh 处理的.」

于是我又去看了与 Dynamic 和 Refresh 有关的内容,这里摘抄其中一条与上文有些关联的例子:

    一个可能会使人伤脑筋的情况是 RandomReal. 每一次你计算 RandomReal[](随机实数 []),你都得到一个不同的答案,你也许就会认为 Dynamic[RandomReal[]](动态[随机实数 []])因此应该不停地尽快地更新自己. 但是在一般情况下这不会很有用,而且实际上会对一些在内部使用随机性的算法有负面的后果(例如,一个在 Dynamic 里面的 Monte Carlo 积分很可能不应该不停地更新,因为事实上它会每次给出一个稍微有些不同的答案).

        Dynamic[Refresh[RandomVariate@NormalDistribution[0, 1], UpdateInterval-> 1]]

这次关于 Mathematica 动态功能的探索就到此为止了,最大的获益就是:以后再遇到什么问题,首先去搜强大的 Help 文档。

2012-07-04

SPICE & LIKES

 

由于原文是 PDF 格式,且有众多公式,不方便转成网页,遂只好以图片形式贴出,哪位朋友如需原文可以把你的邮箱留下,我会将 PDF 文档发给你。

 

 

2011 SPICE & LIKES 综述_页面_1

2011 SPICE & LIKES 综述_页面_2

2011 SPICE & LIKES 综述_页面_3

2011 SPICE & LIKES 综述_页面_4

2012-07-03

Sparse Covariance Fitting for DoA

 

由于原文是 PDF 格式,且有众多公式,不方便转成网页,遂只好以图片形式贴出,哪位朋友如需原文可以把你的邮箱留下,我会将 PDF 文档发给你。

 

 

2012 Covariance Fitting for DoA_页面_1

2012 Covariance Fitting for DoA_页面_2

2012 Covariance Fitting for DoA_页面_3

2012 Covariance Fitting for DoA_页面_4

2012 Covariance Fitting for DoA_页面_5

2012-06-05

Covariance matching estimation techniques for array signal processing applications

由于原文是 PDF 格式,且有众多公式,不方便转成网页,遂只好以图片形式贴出,哪位朋友如需原文可以把你的邮箱留下,我会将 PDF 文档发给你。

 

《Covariance matching estimation techniques for array signal processing applications》加权协方差拟合准则推导_页面_1

《Covariance matching estimation techniques for array signal processing applications》加权协方差拟合准则推导_页面_2

《Covariance matching estimation techniques for array signal processing applications》加权协方差拟合准则推导_页面_3

《Covariance matching estimation techniques for array signal processing applications》加权协方差拟合准则推导_页面_4

《Covariance matching estimation techniques for array signal processing applications》加权协方差拟合准则推导_页面_5

《Covariance matching estimation techniques for array signal processing applications》加权协方差拟合准则推导_页面_6

2012-05-14

LyX 备忘录 - 1

LyX是个很便捷的学术文档编辑器,以下是我最近使用的一些心得(以后会定期更新):

 

1.输入公式可以从 MathType 直接粘贴 latex 代码过来  (最好采用 AMSLaTeX 格式——在 MathType 的粘贴选项可以调整, 并且采用 inline 形式)
然后只需选中代码,按 Ctrl+M 快捷键即可。

 

2.多行公式采用 align 环境对其等号, 逐行粘贴;
把等号以及等号右边的内容放入第二列,第一列留空;
每行还可以单独编号,另外可以分别插入标签、进行交叉引用,很方便。

 

3.按章节编号仅需在导言内写入
\numberwithin{equation}{section}%设置公式按章节进行编号

生成的 latex 代码可以通过 latex2wp.py 程序转换为 Wordpress.com 识别的 HTML (不过正确率不高),方法是:
打开 Windows 命令窗,输入
python “C:\Python26\latex2wp.py" "***\***.tex"

 

2012-05-09

Mathematica Tips – 2

1.Module[]中注意加分号, 否则在计算时会出现意想不到的错误;而且格式化代码时会乱掉.


2.提取列表中的某一项元素出现的位置
Position[{a, b, a, a, b, c, b}, b]


3.Reap 和 Sow 用法:
In[3]:= Reap[Sow[1, {x, x}]; Sow[2, y]; Sow[4, x], {x, x, y}]
Out[3]= {4, {{{1, 1, 4}}, {{1, 1, 4}}, {{2}}}}

关于Reap和Sow以及Scan的一个巧妙用法:
In[18]:= partition[l_, v_, comp_] :=
Flatten /@
  Reap[Scan[
     Which[comp[#1, v], Sow[#1, less], comp[v, #1], Sow[#1, large],
       True, Sow[#1, equ]] &, l], {large, equ, less}][[2]]

In[19]:= partition[{3, 5, 7, 9, 2, 4, 6, 8, 3, 4}, 4, Less]
Out[19]= {{5,7,9,6,8},{4,4},{3,2,3}}


4.不要轻易加脚标:
Subscript[s, k] = Array[0 &, nK];
Subscript[s, k][[1]] = 1;


5.Map(/@) 和 Apply(@@) 对比:

Map不是替换,而是把函数映射到列表的每一个元素, 即把元素分别传递给函数的变量:
In[295]:= f /@ {1, 2, 3, 4, 5}
Out[295]= {f[1], f[2], f[3], f[4], f[5]}

Apply 自动把列表中的元素传递给多变量函数(也就是它仅仅替换函数头),@@@形式自动作用于第一层
In[52]:= Mod[#1, #2] & @@@ {{10, 4}, {5, 2}}
Out[52]= {2,1}

同样的效果, 二者结合方式(借助纯函数):
In[41]:= Apply[Mod, #] & /@ {{10, 4}, {5, 2}}
Out[41]= {2,1}

6.Evaluate函数很有用:
    In[1]:= ch = ChebyshevT[5, x]
    Out[1]= 16 x^5-20 x^3+5 x

    In[4]:= Function[x, Evaluate[ch]]
    Out[4]= x\[Function]16 x^5-20 x^3+5 x

    In[5]:= %[10]
    Out[5]= 1580050
   
   
7.Compile 函数可以提高数值运算速度! (当然采用内部函数是最快的)
它不但可以处理数学表达式,还可以处理各种简单的 Mathematica 程序. 例如,Compile 可以处理条件和控制流结构.
对比:
1.使用Compile:
In[40]:= newtonIteration :=
  Compile[{x, {n, _Integer}}, Module[{t}, t = x;
    Do[t = (t + x/t)/2, {n}]; t]];
newtonIteration[2.4, 666666] // Timing

Out[41]= {0.094,1.54919}

2.不使用Compile:
In[42]:= newtonIteration2[x_, n_] := Module[{t}, t = x;
   Do[t = (t + x/t)/2, {n}]; t];
newtonIteration2[2.4, 666666] // Timing

Out[43]= {2.137,1.54919}

运行时间相差20倍.


8.编写高效代码最重要的方式之一就是要避免显式部分引用,特别是在内部循环中.

如果将要进行实数的操作,一定要确保使用实数进行初始化.

混合的符号/数值矩阵通常比操作数值矩阵要慢.

一个整数矩阵将使用符号计算技术,其速度较慢, 但可以给出精确解。


9.条件迭代:
ff = # + 1 &;
NestWhileList[ff, 2, # < 5 &]

双层迭代的一个例子:
Mp = Table[1, {1}, {4}]
Module[{i, k},
For[i = 1, i < 5, i++,
   Mpp = {};
   Do[AppendTo[Mpp, Sin[k + i]];
    , {k, 4}];
   AppendTo[Mp, Mpp]];]
Mp

 

10.迭代例子2:

intN = 10;
intK = 200;
(*注意是矩阵*)pki = RandomInteger[{1, 2 intK}, {1, intK}];
mA = RandomInteger[{1, 2 intK}, {intK, intN}];
ak = RandomInteger[{1, 2 intK}, {intK, intN}];
y = RandomInteger[{1, 2 intK}, intN];
(*mP[list_]:=DiagonalMatrix[list];*)
(*Subscript[R, y](i)更新*)
Ryi[i_] := mA\[HermitianConjugate].(DiagonalMatrix[pki[[i]]]).mA;
Ryi[1] = RandomInteger[{1, 2 intN}, {intN, intN}];
(*Subscript[p, k](i)更新*)
pkki[k_, i_] := (Abs[ak[[k]]\[Conjugate].Inverse[Ryi[i - 1]].y])^2/(ak[[k]]\[Conjugate].Inverse[Ryi[i - 1]].ak[[k]])^2;


For[i = 1; endCond = 1, endCond > 10^-4(*i<=11*), i++,
bufMp = {};
Do[AppendTo[bufMp,(*N[ pkki[k,i+1] ]*)pkki[k, i + 1]];
  , {k, intK}];
(*AppendTo使用需要注意层次*)
AppendTo[pki, bufMp];
endCond = Norm[pki[[i + 1]] - pki[[i]]]/Norm[pki[[i]]];
]
vectorPk = pki[[i]]

2012-03-22

Mathematica Tips - 1

    注意空格,养成书写严格形式的习惯!
   
    (*不要轻易使用MatrixForm命令*)
    (*使用Take需注意层数*)
   
    局部变量(染成绿色):
        In[2]:= m = i^2
        Out[2]= i^2

        在块内 i+m 的计算过程中,i 用了局部值.
        In[3]:= Block[{i = a}, i + m]
        Out[3]= a + a^2

        这里仅明显出现在 i+m 中的 i 被当作局部变量处理.
        In[4]:= Module[{i = a}, i + m]
        Out[4]= a + i^2

       
    ->和:>区别:
        -> 在首次输入时计算;:> 在使用时计算:
        In[6]:= {x, x, x, x} /. x -> RandomReal[]
        Out[6]= {0.339388, 0.339388, 0.339388, 0.339388}

        In[5]:= {x, x, x, x} /. x :> RandomReal[]
        Out[5]= {0.715386, 0.330938, 0.828099, 0.762862}

   
    Mathematica将向量、矩阵均视为张量处理;与MatLab不同。
       
    数值计算的结果精度不比在变量精度的基础上判断的结果高!
   
    随机数设定种子!
   
   
    @@用来替换函数头部:
        In[2]:= f @@ {a, b, c, d}
        Out[2]= f[a, b, c, d]
       
    /@用来映射函数体
        In[1]:= Reverse /@ {{a, b}, {c, d}, {e, f}}
        Out[1]= {{b, a}, {d, c}, {f, e}}
   
    @,@@,/@区别:
        In[4]:= f@g[a, b, c, d]
        Out[4]= f[g[a, b, c, d]]
       
        In[5]:= f @@ g[a, b, c, d]
        Out[5]= f[a, b, c, d]
       
        In[6]:= f /@ g[a, b, c, d]
        Out[6]= g[f[a], f[b], f[c], f[d]]
       
   
    /.用来替换表达式中的某个变量
        In[2]:= {x, x^2, y, z} /. x -> {a, b}
        Out[2]= {{a, b}, {a^2, b^2}, y, z}
   
    f[##]&    纯函数,##:所有参量
        In[5]:= f[##] & /@ {a, b, c}等价于In[5]:= f[##] & [Sequence[{a, b, c}]]
        Out[5]= {f[a], f[b], f[c]}

        In[4]:= f[##] &[a, b, c]
        Out[4]= f[a, b, c]
   
    MapIndexed传递2个参数:[list,{list num}]
   
    表达式中可以用Print 函数输出:
    Do[If[PrimeQ[2^n - 1], Print[n]], {n, 50, 100}]
   
    Show可以显示多个图形
        Show[g1,g2,...]
       
    Nest可以重复调用函数:
        Nest[Framed, x, 6]
       
    [[]]与[[{}]]的区别,后者有「头部」
    要特别注意「层数」的概念!

2011-05-07

"People Hearing Without Listening:" An Introduction To Compressive Sampling「听而不闻」——压缩采样介绍

本文Google Docs地址

ResearchBlogging.org
Candes, E., & Wakin, M. (2008). An Introduction To Compressive Sampling IEEE Signal Processing Magazine, 25 (2), 21-30 DOI: 10.1109/MSP.2007.914731


1. 介绍


传统信号或者图像采样多采取香农采样定理:即采样频率为信号频率最大值(那奎斯特采样频率)的二倍。

本文提出一种新的信号处理思路——压缩采样(CS)。在这种新方法下,采样率会大大降低。

压缩采样依赖两条原则:稀疏性和不相干性,前者从属于信号,后者从属于感测模式。

  • 稀疏性(sparsity):表达的概念是连续时间信号的信息采集率可能要比按带宽选择的采集率小得多,或者说离散时间信号依赖的自由度数要比他的长度小的多,更明确的说,压缩采样探索如下事实:许多自然信号在下述意义下是稀疏的或可压缩的:当用适当的基表示时,它们有更简洁的表达。

  • 不相干性(incoherence):延伸了时间与频率之间的对偶性,它表达的概念是:正如时间域内的Dirac和尖峰信号在频率域内是展开的那样,具有以$latex {\Psi}&fg=000000$稀疏表示的对象,在获取它们的区域内一定是展开的。


极为重要的是人们可以设计出有效的传感或采样方案,捕捉嵌入在稀疏信号内有用的信息内容,并将其压缩成少量数据。这些方法不是自适应的,且只需要与少量的与提供信号简洁表述之基不相干(否则将失去测量过程对信号的依赖性)的固定波形相关联;此外,有利用数值优化由采集的少量数据重构全长信号的方法。换言之,CS是一种非常简单有效的信息捕获方法,它以与信号无关的方式和低采样率采样,而后根据看似不完整的测量数据集用计算能力重构信号。

2. 感知问题

对于信号$latex {f\left(x\right)}&fg=000000$,传感机制为:

$latex \displaystyle y_{k}=\left\langle f,\varphi_{k}\right\rangle ,\qquad k=1,\ldots,m.\ \ \ \ \ (1)&fg=000000$






也就是说我们只需要将要获得的对象与波形$latex {\varphi_{k}}&fg=000000$相关联,这是一个标准架构,例如,如果检测波形是Dirac delta(尖峰)函数,那么$latex {y}&fg=000000$就是$latex {f}&fg=000000$在时间或空间一个采样值的矢量;如果检测波形是像素的指示函数,那么$latex {y}&fg=000000$就是数字相机中传感器特别采集的图像数据;最后,如果检测波形是正弦函数,那么$latex {y}&fg=000000$就是Fourier系数,磁共振成像 (MRI)用的就是这种检测模态。

尽管可以建立起时空连续的CS理论,但这里只关心离散信号$latex {f\in\mathbb{R}^{n}}&fg=000000$,原因有二,一是简单,二是现有理论在此情况下成熟得多。我们关心欠采样的情况,可得到的测量数$latex {m}&fg=000000$远小于信号$latex {f}&fg=000000$的维数$latex {n}&fg=000000$。由于各种原因,这种情形极为常见,譬如传感器数量有限,或者某些借助于中子散射获取图像的花费极为昂贵,或者是检测过程缓慢,象MRI那样只能对对象施行几次测量,等等。

这些情况产生一些重要问题,根据$latex {m\ll n}&fg=000000$次测量能够准确的重构信号吗?能够设计出$latex {m\ll n}&fg=000000$个检测波形捕获到几乎所有的关于$latex {f}&fg=000000$的信息吗?以及如何根据这些信息逼近$latex {f}&fg=000000$?毫无疑问,这些都是困难的事情,因为可能需要解欠定线性方程组。令$latex {A}&fg=000000$是以$latex {\varphi_{1}^{*},\varphi_{2}^{*},\varphi_{3}^{*},...,\varphi_{m}^{*}}&fg=000000$为行的$latex {m\times n}&fg=000000$测量矩阵,$latex {\varphi^{*}}&fg=000000$是$latex {\varphi}&fg=000000$的复转置。当$latex {m<n}&fg=000000$时,由$latex {y=Af\in\mathbb{R}^{m}}&fg=000000$复原$latex {f\in\mathbb{R}^{n}}&fg=000000$的过程一般不是适定的,有无穷多组解,也许人们会根据具体问题想出一种好办法,这就是本文要讨论的问题。

3. 稀疏性和不相干性

3.1. 稀疏性

通过选择合适的基,许多自然信号都有简洁表示,考虑如图1的图像及其小波变换,尽管原始图像几乎所有的像素都是非零的,但小波系数提供的简明汇总是:绝大部分小波系数的值是小的,为数不多的大系数捕捉了有关对象的主要信息。

用数学的语言描述:我们有一个矢量$latex {f\in\mathbb{R}^{n}}&fg=000000$,以正交基(譬如小波基)$latex {\Psi=\left[\Psi_{1},\Psi_{1},\ldots,\Psi_{1}\right]}&fg=000000$展开如下:

$latex \displaystyle \ensuremath{{\rm {f}}\left({\rm {t}}\right)=\mathop\sum\limits _{{\rm {i}}=1}^{{\rm {n}}}{{\rm {x}}_{{\rm {i}}}}{{\rm {\psi}}_{{\rm {i}}}}({\rm {t}})}\ \ \ \ \ (2)&fg=000000$






$latex {x}&fg=000000$是系数列,$latex {x_{i}=\left\langle f,\Psi_{i}\left(t\right)\right\rangle }&fg=000000$,将$latex {f}&fg=000000$表示成$latex {\mathbf{\Psi x}}&fg=000000$。现在明确稀疏性的含义:当信号有稀疏展开时,可以丢掉小系数而不会失真。按正式的说法,$latex {f_{s}(t)}&fg=000000$是保留展开式(2)中$latex {S}&fg=000000$个最大系数$latex {(x_{i})}&fg=000000$值得到的结果。有定义$latex {f_{s}:=\Psi x_{s}}&fg=000000$,此后$latex {x_{s}}&fg=000000$就是系数向量$latex {x_{i}}&fg=000000$,只不过除了$latex {S}&fg=000000$个最大值外其余都是0,该向量在严格意义上是稀疏的;因为,除了少数几个非零元素外,其余元素都为0。

我们称这种几乎$latex {S}&fg=000000$个非零元的对象为$latex {S}&fg=000000$-稀疏的,由于$latex {\Psi}&fg=000000$是正交基,我们有:$latex {\|f-f_{s}\|_{\ell_{2}}=\|x-x_{s}\|_{\ell_{2}}}&fg=000000$

如果$latex {x}&fg=000000$在按值排序快速衰减的意义上是稀疏的,$latex {x}&fg=000000$就能很好地用$latex {x_{s}}&fg=000000$逼近,误差$latex {\|f-f_{s}\|_{\ell_{2}}}&fg=000000$就是小量。通俗点说就是除了几个大系数外扔掉其余系数,不会造成太大的损失。从图1给出的例子,几乎觉察不到1兆像素图像与丢掉97.5%系数的近似图像之间的差别。


这一原理为现代有损失编码基础,如JPEG-2000及其它编码器,由于有了一种数据压缩的简单方法,将只需根据$latex {f}&fg=000000$计算$latex {x}&fg=000000$,然后(自适应地)将$latex {S}&fg=000000$个重要系数的值及位置进行编码。这种压缩方法需要知道所有$latex {n}&fg=000000$个系数$latex {x}&fg=000000$,由于事先不知道重要信息片段的位置(它们与信号有关,在我们的例子中有积聚在图像边缘的趋势,其位置事先未知)。在更一般的意义上,稀疏性是一种基本模型工具,它允许有效的基本信号处理,例如:准确的统计计算和分类,有效的数据压缩等等。本文要研究的事情是它更惊奇深远的含义,这就是信号稀疏性在数据采集本身所具有的重要作用。稀疏性决定着人们如何能够有效、非自适应的获取信号。

3.2. 非相干采样

假定我们给定$latex {\mathbb{R}^{n}}&fg=000000$内一组正交基$latex {\left(\Phi,\Psi\right)}&fg=000000$,第一个基$latex {\Phi}&fg=000000$用于像式(1)那样感知对象$latex {f}&fg=000000$,第二个基用于表示$latex {f}&fg=000000$,对这一对基的限制不是实质性的,仅仅为了简化处理。

定义1 :

$latex {\Phi}&fg=000000$和$latex {\Psi}&fg=000000$之间的相干度:

$latex \displaystyle \mu\left(\Phi,\Psi\right)=\sqrt{x}\mathop{\max}\limits _{1\le k,j\le n}\left|{\left\langle {{\varphi_{k}},{\psi_{j}}}\right\rangle }\right|\ \ \ \ \ (3)&fg=000000$






用通俗的话说,相干度度量着$latex {\Phi}&fg=000000$和$latex {\Psi}&fg=000000$中任意两个元素的最大相关性。如果$latex {\Phi}&fg=000000$和$latex {\Psi}&fg=000000$含有相干元,相干系数$latex {\mu}&fg=000000$值就很大,否则值是小量。无论大和小要满足$latex {\mu\left(\Phi,\Psi\right)\in\left[1,\sqrt{n}\right]}&fg=000000$。上端点是因为基为单位向量,下端点是由于Parseval关系的结果,该关系指出,对每一个$latex {j}&fg=000000$,$latex {\sum\nolimits _{k=1}^{n}{{{\left|{\left\langle {{\varphi_{k}},{\Psi_{j}}}\right\rangle }\right|}^{2}}=\left\Vert {\Psi_{j}}\right\Vert _{{\ell_{2}}}^{2}}=1}&fg=000000$。压缩采样主要与低相干度对有关,下面给出这样相干度对的一些例子。

第一个例子中$latex {\Phi}&fg=000000$是经典的或尖峰基$latex {\varphi_{k}\left(t\right)=\delta\left(t-k\right)}&fg=000000$,$latex {\Psi}&fg=000000$是Fourier基,$latex {{\Psi_{j}}\left(t\right)={n^{{{-1}\mathord{\left/{\vphantom{{-1}2}}\right.\kern -\nulldelimiterspace}2}}}{e^{{{i2\pi jt}\mathord{\left/{\vphantom{{i2\pi jt}n}}\right.\kern -\nulldelimiterspace}n}}}}&fg=000000$。由于$latex {\Phi}&fg=000000$是测量矩阵,这相应于时、空经典的采样格式,时间频率对遵循$latex {\mu\left(\Phi,\Psi\right)=1}&fg=000000$,因此有最大相干性。此外,尖峰信号与正弦波不仅在1-D情况下不相干,在任意维数都如此,2D、3D等等相同。

第二个例子$latex {\Psi}&fg=000000$是小波基,$latex {\Phi}&fg=000000$取Noiselets。Noiselets同Harr小波基的相干度是$latex {\sqrt{2}}&fg=000000$,同Daubechies D4和 D8小波的相干度分别是2.2和2.9,跨越非常大的样本容量$latex {n}&fg=000000$的范围;延伸至高维情况也是如此。(Noiselets同尖峰信号也是最大不相干的,同Fourier基不相干)。由于以下事实我们对Noiselets感兴趣:(1)它们与提供图像数据及其它数据稀疏表示的系统不相干;(2)它们适合于非常快速的算法,Noiselets变换耗时$latex {O(n)}&fg=000000$, 和Fourier变换一样Noiselets矩阵再作用于向量时不需存储,这对于有效的数值计算是极为重要的,没有这一点CS不可能很实用。

最后一个例子,随机矩阵基本上同任何固定的基$latex {\Psi}&fg=000000$都是不相干的。可以通过在单位球上独立均匀采样的$latex {n}&fg=000000$个向量正交化,均匀随机地 选一个正交基$latex {\Phi}&fg=000000$,那么,以高概率使$latex {\Phi}&fg=000000$和$latex {\Psi}&fg=000000$相干度大约为$latex {\sqrt{2\log n}}&fg=000000$。通过扩充具有独立同分布元素的随机波形$latex {\varphi_{k}\left(t\right)}&fg=000000$也将展示出与固定表示$latex {\Psi}&fg=000000$有非常低的相干性,例如Gauss型或$latex {\pm1}&fg=000000$二进制元。注意到这里有一个非常奇特的暗示:如果不相干系统测量是良好的,那么高效的机制应该获得同随机波形的关联,例如白噪声!

3.3. 欠采样和稀疏信号重建

从理论上,我们本来测的是$latex {f}&fg=000000$的$latex {n}&fg=000000$个系数,但是我们只观测到它们的子集,并采集数据:

$latex \displaystyle y_{k}=\left(f,\varphi_{k}\right),\quad k\in M\ \ \ \ \ (4)&fg=000000$






这里$latex {M\subset\left\{ 1,\ldots,n\right\} }&fg=000000$是基数$latex {m<n}&fg=000000$的子集,用这一信息,我们决定由$latex {\ell_{1}}&fg=000000$-范数($latex {\|x\|_{\ell_{1}}:=\sum_{i}|x_{i}|}&fg=000000$)极小化复原信号,所提出的重构$latex {f^{*}}&fg=000000$由$latex {f^{*}=\Psi x^{*}}&fg=000000$给出;$latex {x^{*}}&fg=000000$是下述凸优化问题的解:

$latex \displaystyle {\min_{\tilde{x}\in{\mathbb{R}^{n}}}}{\left\Vert {\tilde{x}}\right\Vert _{{\ell_{1}}}}\qquad s.t.\qquad{y_{k}}=\left\langle {{\varphi_{k}},\Psi\tilde{x}}\right\rangle ,\quad\forall k\in M\ \ \ \ \ (5)&fg=000000$






即,在所有与数据相容的对象$latex {f^{*}=\Psi x^{*}}&fg=000000$中挑选出最小$latex {\ell_{1}}&fg=000000$-范数的系数序列。

用$latex {\ell_{1}}&fg=000000$-范数作为稀疏提升函数可以追溯回几十年前。早先应用在反射地震学,从限带宽数据中寻找稀疏反射函数(用以指示地表下各层的重要变化)。但$latex {\ell_{1}}&fg=000000$-极小化不是复原稀疏解的唯一方法,也提出了其它方法,如贪婪算法。

第一个结果断言:当$latex {f}&fg=000000$足够稀疏时,可以证明,籍$latex {\ell_{1}}&fg=000000$-极小化复原信号是准确的。

定理1 :

固定$latex {f\in\mathbb{R}^{n}}&fg=000000$,并假设基$latex {\Psi}&fg=000000$下系数序列$latex {x}&fg=000000$是$latex {S}&fg=000000$-稀疏的,在$latex {\Phi}&fg=000000$域上随机均匀地选择$latex {m}&fg=000000$次测量。那么,如果对某一正常数$latex {C}&fg=000000$有:

$latex \displaystyle m\geq C\mu^{2}\left(\Phi,\Psi\right)\cdot S\cdot\log n\ \ \ \ \ (6)&fg=000000$






(5)式的准确率是可以充分保证的。

这里做三点解释:

  1. 相干性的作用是清晰的,相干度越小,需要的采的样就越少,因此,在前面一节中强调低相干度;

  2. 通过测量任一组$latex {m}&fg=000000$系数(可能远小于信号的表面需要)人们都不会蒙受信号损失;如果$latex {\mu\left(\Phi,\Psi\right)}&fg=000000$等于或接近于1,那么$latex {S\cdot\log n}&fg=000000$次采样就足够了,而不是$latex {n}&fg=000000$次;

  3. 在事先不知道$latex {x}&fg=000000$非零坐标个数、位置的条件下,信号$latex {f}&fg=000000$可通过凸泛函极小化根据压缩数据集复原,关于它们的幅值事先完全未知。


这个定理确实提出了一个非常具体的捕获方案,在不相干域采样不是自适应的,采集后执行线性规划。按这一方法,得到实质上的压缩信号,所需要的是将数据解压缩的解码器,这是$latex {\ell_{1}}&fg=000000$-极小化所起的作用。

事实上这个非相干采样是早先谱稀疏信号采样结果的推广;也许正是谱稀疏采样引发了我们已经证明而且今天继续证明的CS发展。假设我们对超宽带采样感兴趣,而谱稀疏信号形为:

$latex \displaystyle f\left(t\right)=\sum_{j=0}^{n-1}x_{j}e^{i2\pi jt/n},\, t=0,\dots,n-1,&fg=000000$


此处虽然$latex {n}&fg=000000$很大,但非零分量$latex {x_{j}}&fg=000000$小于或等于$latex {S}&fg=000000$(我们可以理解为非常小)。我们不知道哪些频率是主要的,也不知道主要频率集合上的振幅。由于主分量集不一定是顺序整数的子集,Nyquist/Shannon理论几乎对我们毫无帮助(因为我们不能事先限制带宽,甚至认为所有$latex {n}&fg=000000$倍采样都需要)。就这个特例而言,定理1告诉人们可以根据$latex {S\cdot\log n}&fg=000000$采样顺序,用任意未知容量$latex {S}&fg=000000$的频率支集重构信号;此外,这些采样不必是精心选择的,几乎这个容量里的任一样本集都有效;图2给出了一个说明性的例子。


现在讨论在这方面概率扮演的角色;人们需要借助概率阐述,因为人们不能希望容量$latex {m}&fg=000000$内的所有测量集都成立着类似结果。原因是有一些特殊的稀疏信号几乎在整个$latex {\Phi}&fg=000000$域都为零,换言之,可以找到稀疏信号$latex {f}&fg=000000$和几乎等于$latex {n}&fg=000000$非常大容量的子集(e.g. $latex {n-S}&fg=000000$),对所有$latex {k\text{\ensuremath{\in}}M}&fg=000000$,$latex {y_{k}=\left\langle f,\varphi_{k}\right\rangle =0}&fg=000000$。一方面给出这样的子集人们看到的尽是零,当然没有算法重构这种信号;另一方面,定理保证数据集中不发生准确复原的那一部分实际上是可忽略的;因此我们必须容忍极小的失败概率。对实用而言,如果采样容量足够大,失败的概率是零。

有趣的是,上述讨论的特殊稀疏信号至少也需要$latex {\mu^{2}\cdot S\cdot\log n}&fg=000000$个样本。用小的采样容量只是信息损失的概率太高,用任意多么复杂的方法重构都是不可能的。总的说,当相干度是1时,采样数不必大于$latex {S\cdot\log n}&fg=000000$,但也不能再小。

我们以不相干采样的例子结束这一节。考虑图3所示演员的稀疏图像,与图1压缩的百万像素相同。研究对象是稀疏的,由于它的非零小波系数仅有25000个,然后取96000次非相干测量获取信息,并求解问题(5)。该例表明采样数大约为稀疏水平的4倍就满足了。许多研究者也报道了类似的成功经验,事实上有一个4:1的实用规则,该规则告知,为了准确的复原信号,每个未知非零项需要4个不相关采样。


4. 稳健压缩采样

我们已经证明可以仅由少量测量复原稀疏信号,但为了实际上更为有效起见,CS需要能够处理带有噪声的近稀疏信号。

  1. 首先,一般关心的对象并不是严格而是近似稀疏的,问题是:是否能够从高度欠采样测量中准确的复原这些对象。

  2. 其次,在任何实际应用的测量数据中,由于测量设备不会有无限精度,总是存在哪怕是很小的噪声干扰。最轻微的情况下,噪声对数据的小扰动也会引起对重构的小扰动。


本节同时考察这两个问题。便于叙述,开始之前先考虑从形如

$latex \displaystyle y=Ax+z\ \ \ \ \ (7)&fg=000000$






复原向量$latex {x\in\mathbb{R}^{n}}&fg=000000$,$latex {A}&fg=000000$是告诉我们$latex {x}&fg=000000$信息的$latex {m\times n}&fg=000000$测量矩阵,$latex {z}&fg=000000$是随机的或确定性的误差项。由于$latex {f=\Psi x}&fg=000000$,$latex {y=R\Phi f}&fg=000000$($latex {R}&fg=000000$是从$latex {M}&fg=000000$中提取采样坐标的$latex {m\times n}&fg=000000$矩阵),令$latex {A=R\Phi\Psi}&fg=000000$,可以写成$latex {y=Ax}&fg=000000$,因此人们可以用抽象模型(7)。$latex {x}&fg=000000$可以是以适当的基表示对象的系数序列。

4.1. 约束等距性

约束等距性:restricted isometry property (RIP),是对CS的一般鲁棒性非常有用的一个重要概念。

定义2:

对每一个整数$latex {s=1,2,\ldots,}&fg=000000$定义矩阵A的等距常数$latex {\delta_{s}}&fg=000000$为使下式:

$latex \displaystyle \left({1-{\delta_{s}}}\right)\left\Vert x\right\Vert _{{\ell_{2}}}^{2}\leqslant\left\Vert {Ax}\right\Vert _{{\ell_{2}}}^{2}\leqslant\left({1+{\delta_{s}}}\right)\left\Vert x\right\Vert _{{\ell_{2}}}^{2}\ \ \ \ \ (8)&fg=000000$




对所有$latex {s}&fg=000000$-稀疏向量$latex {x}&fg=000000$成立的最小数。

如果$latex {\delta_{s}}&fg=000000$不太接近于1,就不十分严谨的说矩阵$latex {A}&fg=000000$遵循$latex {s}&fg=000000$阶约束等距性。如果这一性质成立,$latex {A}&fg=000000$近似保持$latex {S}&fg=000000$-稀疏信号的欧几里得长度,反过来说$latex {S}&fg=000000$-稀疏向量不可能在$latex {A}&fg=000000$的零空间(这是有用的,若不其然就没有重构信号的希望)。约束等距性的等价表述是,取自$latex {A}&fg=000000$的$latex {S}&fg=000000$列子集事实上近乎正交($latex {A}&fg=000000$的列不可能准确正交,因为行大于列)。

为看到约束等距与CS之间的联系,想象一下我们用$latex {A}&fg=000000$获取$latex {S}&fg=000000$-稀疏信号。首先假设$latex {\delta_{2S}<1}&fg=000000$,我们就声称能由数据$latex {y=Ax}&fg=000000$复原信号$latex {x}&fg=000000$,事实上$latex {x}&fg=000000$是方程组$latex {y=A\tilde{x}}&fg=000000$的唯一稀疏解,即非零元的个数最少。为看清这一点,考虑任意其它解$latex {x+h,h\neq0}&fg=000000$,那么$latex {Ah=0}&fg=000000$,$latex {h}&fg=000000$一定至少有$latex {2S+1}&fg=000000$个非零元素,则$latex {x+h}&fg=000000$至少有$latex {S+1}&fg=000000$个非零元。相反,假设$latex {\delta_{2S}=1}&fg=000000$,那么,$latex {A}&fg=000000$的2S列必定线性相关,在此情形下2S-稀疏向量$latex {h}&fg=000000$遵循$latex {Ah=0}&fg=000000$,可以将$latex {h}&fg=000000$分解为$latex {x-x'}&fg=000000$使$latex {x}&fg=000000$和 $latex {x'}&fg=000000$都是$latex {S}&fg=000000$-稀疏信号,由此得出$latex {Ax-Ax'}&fg=000000$,这意味着我们找到了给出相同测量的两个稀疏向量,显然人们不能重构这样的稀疏信号。由此可见,要复原$latex {S}&fg=000000$-稀疏信号,必须满足$latex {\delta_{2S}<1}&fg=000000$,或至少(8)中$latex {\|Ax\|_{\ell_{2}}}&fg=000000$取低限。

4.2. 欠采样复原信号

如果满足约束等距性,解线性规划问题:

$latex \displaystyle {\min_{\tilde{x}\in{\mathbb{R}^{n}}}}{\left\Vert {\tilde{x}}\right\Vert _{{\ell_{1}}}}\qquad s.t.\qquad A\tilde{x}=y\left({=Ax}\right)\ \ \ \ \ (9)&fg=000000$




得到的重构是准确的。

定理 2

假设$latex {\delta_{2S}<\sqrt{2}-1}&fg=000000$,那么问题(9)的解$latex {x^{*}}&fg=000000$对某一常数$latex {C_{0}}&fg=000000$遵循:

$latex \displaystyle {\left\Vert {{x^{*}}-x}\right\Vert _{{\ell_{2}}}}\leqslant{C_{0}}\cdot{\left\Vert {x-{x_{S}}}\right\Vert _{{\ell_{1}}}}/\sqrt{S}\quad and\quad{\left\Vert {{x^{*}}-x}\right\Vert _{{\ell_{1}}}}\leqslant{C_{0}}\cdot{\left\Vert {x-{x_{S}}}\right\Vert _{{\ell_{1}}}}\ \ \ \ \ (10)&fg=000000$




$latex {x_{S}}&fg=000000$是除去最大的$latex {S}&fg=000000$个值其余置为零分量的向量$latex {x}&fg=000000$。这个定理的结论比定理1强。若$latex {x}&fg=000000$是$latex {S}&fg=000000$-稀疏的,$latex {x=x_{S}}&fg=000000$,信号复原是准确的;若$latex {x}&fg=000000$不是$latex {S}&fg=000000$-稀疏的,(10)式断言:复原信号的质量恰如人们提前知道最大$latex {S}&fg=000000$个最大$latex {x}&fg=000000$值并确定直接测量的结果那么好。换言之,重构近乎对对象完全了解的预言,为我们抽取出$latex {S}&fg=000000$个最重要的信息片段。

与我们早先的结果的另一个惊人的差别在于它是确定性的,不涉及到概率。如果我们足够幸运使测量矩阵$latex {A}&fg=000000$遵循定理的假设,就可以用它,并保证准确地复原所有$latex {S}&fg=000000$-稀疏信号,而对其它情况实质上是复原所有矢量的$latex {S}&fg=000000$个最大元素;即,没有失败的概率。此外,相同测量矩阵对$latex {\mathbb{R}^{n}}&fg=000000$中所有向量都适用,常称为通用方法。

用此方法我们失去的是遵循假设的$latex {S}&fg=000000$(能有效复原的分量数)与$latex {m}&fg=000000$(测量数或矩阵的行)之间的关系,为了导出更有效的结果,我们宁愿寻找$latex {S}&fg=000000$接近$latex {m}&fg=000000$满足约束等距性的测量矩阵。能设计出这样的矩阵吗?我们下一节将证明这是可能的,但首先考察一下面对数据干扰CS的鲁棒性。

4.3. 噪声信号的稳健恢复

给我们如(7)所示的有噪声数据,用$latex {\ell_{1}}&fg=000000$极小化重构约束:

$latex \displaystyle \min{\left\Vert {\tilde{x}}\right\Vert _{{\ell_{1}}}}\qquad s.t.\qquad{\left\Vert {A\tilde{x}-y}\right\Vert _{{\ell_{2}}}}\leqslant\epsilon\ \ \ \ \ (11)&fg=000000$






$latex {\epsilon}&fg=000000$为数据噪声值的界,问题(11)在文献[21]、[22]中称为LASSO;就我们所知这首先是在[8]中提出的。这又是一个特有的凸优化问题,准确的说是二阶圆锥规划,求解此问题有多种有效的算法。

定理 3:

假设$latex {\delta_{2S}<\sqrt{2}-1}&fg=000000$,那么问题(11)的解$latex {x^{*}}&fg=000000$对某一常数$latex {C_{0}}&fg=000000$和$latex {C_{1}}&fg=000000$遵循:

$latex \displaystyle {\left\Vert {{x^{*}}-x}\right\Vert _{{\ell_{2}}}}\leqslant{C_{0}}\cdot{\left\Vert {x-{x_{S}}}\right\Vert _{{\ell_{1}}}}/\sqrt{S}+{C_{1}}\cdot\epsilon\ \ \ \ \ (12)&fg=000000$






这几乎不能更简单了,重构误差有两项之和限定,第一项是没有噪声时的误差,第二项正比于噪声水平。此外,$latex {C_{0}}&fg=000000$和$latex {C_{1}}&fg=000000$是特征小量,例如在$latex {\delta_{2S}=1/4}&fg=000000$的情况下,$latex {C_{0}\leq5.5}&fg=000000$和$latex {C_{1}\leq6}&fg=000000$。图4显示有噪声数据的重构。


最后这个结果为CS建立了实用有效的传感机制。它对各种不一定是稀疏的信号都有效,在传感方面能优雅地处理噪声。剩下要做的事情就是设计满足约束等距性(RIP)的测量矩阵;这是下一节的目标。

5. 随机测量

回到RIP的定义,我们要寻找测量矩阵具有列矢量的任意子集都几乎正交的性质;这个子集越大越好;为此,随机性重新出现。

考虑如下测量矩阵:

  1. 通过在单位球$latex {\mathbb{R}^{m}}&fg=000000$上随机均匀地采集$latex {n}&fg=000000$列矢量,形成$latex {A}&fg=000000$。

  2. 通过从平均值为0,方差为$latex {1/m}&fg=000000$的正态分布上独立同分布(i.i.d.)采集元素,形成$latex {A}&fg=000000$。

  3. 如节3.2所述对随机投影$latex {P}&fg=000000$采样并标准化为:$latex {A=\sqrt{\frac{n}{m}}P}&fg=000000$,形成$latex {A}&fg=000000$。

  4. 从对称Bernoulli分布$latex {\left(P\left(A_{i,j}=\pm1/\sqrt{m}\right)=\frac{1}{2}\right)}&fg=000000$或其它亚Gauss分布上独立同分布采集元素,形成$latex {A}&fg=000000$。


假如:

$latex \displaystyle m\geq C\cdot S\log\left(n/S\right)\ \ \ \ \ (13)&fg=000000$






那么,所有这些矩阵将以极大概率符合约束等距性(即定理的条件),上式中$latex {C}&fg=000000$为与具体情况有关的某一常数。命题1~3用了概率论中非常标准的结果。对命题4的论证更复杂些,见[23]和Pajor及其合作者的研究,如[24]。

在所有这些例子中,测量矩阵满足(13)式而又不具有约束等距性的概率是关于$latex {m}&fg=000000$指数型小量。有趣的是,没有测量矩阵和重构算法用比式(13)左端少得多的采样而又能给出定理2的结果。在此意义上,用上述测量矩阵及$latex {\ell_{1}}&fg=000000$极小化是接近最优的感知策略。

可以像第三节中所作的那样建立正交基对的约束等距性,用$latex {A=R\Phi\Psi}&fg=000000$,$latex {R}&fg=000000$随机均匀地抽取$latex {m}&fg=000000$个坐标,使

$latex \displaystyle m\geq C\cdot S\left(\log n\right)^{4}\ \ \ \ \ (14)&fg=000000$






是以极大概率具有该性质的充分条件。如果要失败的概率对某$latex {\beta>0}&fg=000000$不大于$latex {O\left(n^{-\beta}\right)}&fg=000000$,知(14)中的指数是5而不是4(人们相信(14) 对恰为$latex {\text{\ensuremath{\log}}n}&fg=000000$也是成立的)。这证明可以稳定而准确地由不相干域中戏剧性的欠采样数据重构近稀疏信号。

最后,矩阵$latex {A=\Phi\Psi}&fg=000000$也可以成立约束等距性,这里$latex {\Psi}&fg=000000$是任意正交基,$latex {\Phi}&fg=000000$是从适当分布中随机抽取的测量矩阵。如果固定基$latex {\Psi}&fg=000000$,按命题1~4构成测量矩阵,假若:

$latex \displaystyle m\geq C\cdot\log\left(n/S\right)\ \ \ \ \ (15)&fg=000000$






那么,$latex {A=\Phi\Psi}&fg=000000$将以极大的概率符合约束等距性(RIP);式中$latex {C}&fg=000000$为与具体情况有关的某一常数。这种随机测量矩阵$latex {\Phi}&fg=000000$在某种意义上是普适的,当设计测量系统时甚至不需知道稀疏基。

6. 什么是压缩采样?

数据获取通常按下述方式工作:大量的数据只是采集出来,其中一大部分在压缩阶段被抛弃,为了储存和传送通常必须这样做。用本文的话说,人们获得高分辨率的像素数组$latex {f}&fg=000000$,计算完整的一套变换系数,将最大的系数编码,丢掉其它系数,实质上以$latex {f_{S}}&fg=000000$结束;这种大规模采集数据然后压缩的处理方法是极为浪费的(可以想象一下数字相机有百万像素的图像传感器,而最终编码的图像只有几百k字节)。

CS的处理极为不同,其表现好像是能够直接获取处理对象恰为最重要的信息。通过取如第5节的$latex {O(S\log(n/S))}&fg=000000$个随机投影,人们就有足够的信息重构信号,其精度至少达到$latex {f_{S}}&fg=000000$具有的精度,得到处理对象的最佳$latex {S}&fg=000000$项逼近和最优压缩表示。换言之,CS数据获取协议实质上是把模拟信号转换成压缩后的数字形式,使得人们至少从原理上可以从仅有为数不多的传感器得到超分辨信号。采集步骤后所有需要做的就是将测量数据解压缩。料想不到的是采集步是固定的,尤其是根本不试图理解信号的结构,好像是「听而不闻」( hearing without listening)。

CS同编码理论,更明确的说同Reed-Solomon编码的理论和实践有某些表面上的相似性。就本文讨论内容简单地说,如所周知,可以采用编码理论的概念如下建立CS:人们可以根据信号的前2S个 Fourier系数

$latex \displaystyle {y_{k}}=\sum\limits _{t=0}^{n-1}{{x_{t}}{e^{-i2\pi kt/n}}},\quad k=0,1,2,\ldots,2S-1\ \ \ \ \ (16)&fg=000000$






或者从相继2S 个信号频率集唯一的重构任何$latex {S}&fg=000000$-稀疏信号(复原信号的计算成本实质上是解$latex {S\times S}&fg=000000$ Toeplitz矩阵,或计算$latex {n}&fg=000000$点FFT),这意味着可以用这一方法测量可压缩信号吗?由于两个原因答案是否定的。

其一,Reed-Solomon解码是代数方法,不能用于非稀疏信号(解码是通过求多项式的根寻找信号支撑集);

其二,根据信号的前2S个Fourier系数寻找信号支撑集的问题,甚至当信号具备准确稀疏性时,是非常不适定的(这个问题与由少量值高度集中的数预测高阶多项式相同);哪怕系数受到微小的扰动也会导致完全不同的答案,这使得用有限精度的数据可靠地计算支撑集在实际上是不可能的。反之,纯粹代数方法忽略了信息算子的调节作用,使取得良态矩阵,这是精确估计的关键,正如约束等距性所起作用是CS所关心的中心问题。

7. 应用

可压缩信号可以用与信息水平$latex {S\ll n}&fg=000000$成正比的若干不相干测量捕获,这一事实有着深远的意义并涉及到许多可能的应用。

  1. 数据压缩:对数据压缩而言,在某些情况下解码器上稀疏基可能是未知的,或者施行起来不实用。然而正如第五节所讨论的随机设计的$latex {\Phi}&fg=000000$可视为普适解码策略,因为它不需要围绕$latex {\Psi}&fg=000000$的结构设计(只有解码或复原$latex {f}&fg=000000$时才需具备$latex {\Psi}&fg=000000$的知识和实现$latex {\Psi}&fg=000000$的能力)。这种普适性特别有助于诸如传感器网络之类的多信号装置的编码。关于这个问题请读者参考Nowak等人 及Goyal在其它地方的论文。

  2. 信道编码:正如文献[15]所解释的那样,可以围绕CS的原理(稀疏性,随机性,和凸优化)并将其用于设计快速纠错码,避免错误信号的传输。

  3. 反问题:对于其它一些情况,获取$latex {f}&fg=000000$的唯一方法可能是用某种模态的测量系统$latex {\Phi}&fg=000000$,假定与$latex {\Phi}&fg=000000$不相干$latex {f}&fg=000000$的稀疏基$latex {\Psi}&fg=000000$存在。就有可能进行有效的测量,一个应用涉及到MR血管造影术和其它类型的MR设备;在这些例子中$latex {\Phi}&fg=000000$记录$latex {f}&fg=000000$的Fourier变换的子集,所希望的图像$latex {f}&fg=000000$在时间域和小波域都是稀疏的。其它领域的这个问题Lustig等人有更深入的讨论。

  4. 数据采集:有时全部测量模拟信号的n个离散时间样本可能难以得到(而且接下来也难以压缩)。在这种情况下,可能有助于设计出物理采样设备,直接记录传入模拟信号离散的、低采样率的不相干测量。


最后一个应用暗示我们,数学和计算方法可能会对传统的硬件设计有限制的领域产生巨大的冲击。如传统的图像设备采用CCD和CMOS技术基本上限于可视光谱,而CS相机用数字微型镜头阵列采集不相干数据,也许能够明显扩充其能力。关于这个问题其它地方的讨论,以Baraniuk讨论这种相机更加深入。

有一部分研究者已专心于大带宽高级模拟-信息(A/I)转换设备的研究,目标是帮助减轻传统模拟-数字(A/D)转换技术的压力,当前该技术限于1GHz的采样速度。作为可选择的方案,我们建议两种特殊的、可以从高带宽模拟信号中获得离散的、低采样率,不相干测量序列的A/I架构。对高阶逼近,每一个测量$latex {y_{k}}&fg=000000$,都可以解释为入射模拟信号$latex {f}&fg=000000$与模拟测量波形$latex {\varphi_{k}}&fg=000000$的内积$latex {\left\langle f,\varphi_{k}\right\rangle }&fg=000000$,正如在离散CS框架那样,我们的初步结果表明遵循稀疏或可压缩模型(在某一模拟字典$latex {\Psi}&fg=000000$中)的信号可以用这些设备以正比于信息水平而不是那奎斯特频率有效的捕获。当然,利用离散CS方法复原模拟稀疏信号有待克服的挑战性问题。全面解决这些问题超出了本短文的范围;人们可以简单的接受其思想,离散/采样稀疏字典允许适当的复原。

  • 我们有两个架构如下:



  1. 非均匀采样器(NUS):此架构只是把随机采样时间点上的信号数字化,即$latex {{y_{k}}=f\left({t_{k}}\right)=\left\langle {f,{\delta_{{t_{k}}}}}\right\rangle }&fg=000000$,事实上这些随机或伪随机时间点是通过抖动规则格子上的名义(低采样率)采样点得到的。由于尖峰与正弦信号不相干,因此这种架构可用于对具有远低于那奎斯特频率的稀疏频谱信号采样。当然与减少采样率有关的好处是极大的,因为这提供了附加电路稳定时间,并具有减小噪声水平的作用。

  2. 随机预积分(RPI):此架构可用于种类更广泛的稀疏域,在时间-频率平面有稀疏特征的最显著的那些信号。鉴于没有可能以极高的采样率将模拟信号数字化,而极有可能以很高的速率改变它的极性。RPI架构的思想(见图5)就是以正1和负1的伪随机序列乘以信号,将乘积在时间窗口积分,将时段末的积分值数字化,这是一个并行架构,有若干个这样的乘法器-积分器对用完全不同的甚至几乎独立的随机符号序列并列运行。事实上RPI架构将信号同$latex {\pm1}&fg=000000$符号序列库相关联,是普适的CS测量方法之一。


对上述每个架构,我们都曾用数值的方法确认,有些还用物理的方法证实,系统对于电路非理想,如热噪声、时钟误差、干扰及放大器非线性等情况是鲁棒性的。

将A/I架构应用于实际采集方案还需要发展CS算法和理论,包括研究模拟信号稀疏表示的字典。我们最后以一个离散的例子结束。对这个实验,我们取$latex {f}&fg=000000$是长度$latex {n=512}&fg=000000$的1-D信号,包含两个调制脉冲(图6左侧),从这个信号中我们用以独立同分布柏努利$latex {\pm1}&fg=000000$元素填充生成的$latex {m\times n}&fg=000000$测量矩阵$latex {\Phi}&fg=000000$, 采集$latex {m=30}&fg=000000$次测量,这是不合理的小数据量,欠采样因子超过17。为了重构信号,我们考虑时-频Gabor字典$latex {\Psi}&fg=000000$,它由Gaussian窗限时并具有不同位置和尺度的各种正弦波构成。总的说字典是过完备($latex {43\times}&fg=000000$overcomplete)的,不含包括$latex {f}&fg=000000$的两个脉冲。图6中间图形表示取$latex {{\left\Vert x\right\Vert _{{\ell_{1}}}}}&fg=000000$极小化使$latex {y=\Phi\Psi x}&fg=000000$,重构效果显示出明显的人为效应,我们看到$latex {{\left\Vert {f-{f^{*}}}\right\Vert _{{\ell_{2}}}}/{\left\Vert f\right\Vert _{{\ell_{2}}}}\approx0.67}&fg=000000$,而事实上我们通过将两者都变为$latex {\ell_{1}}&fg=000000$复原方法,就可以消除人为效应。首先,我们改用极小化$latex {{\left\Vert {\Psi*\tilde{f}}\right\Vert _{{\ell_{1}}}}\; s.t.\; y=\Phi\tilde{f}}&fg=000000$;(当$latex {\Psi}&fg=000000$是正交基时这一改变不起作用);其次,当得到估计值$latex {f^{*}}&fg=000000$时我们将$latex {\ell_{1}}&fg=000000$-范数重新加权,重复重构过程,对那些估计较大的系数采用较低的惩罚值;图6右图显示了重复加权4次迭代的结果,我们看到$latex {{\left\Vert {f-{f^{*}}}\right\Vert _{{\ell_{2}}}}/{\left\Vert f\right\Vert _{{\ell_{2}}}}\approx0.022}&fg=000000$。有关这个新方向的更多信息,建议读者参阅文献[30]。在此,给出的要点是,即便这个数据量小到荒谬的地步,人们仍可以捕获到信号里包含的绝大部分信息。这就是CS在当前以及未来应用极有前途的原因。


 

References
[1] E.J. Cand`es, J. Romberg, and T. Tao. Robust uncertainty principles: Exact signal reconstructionfrom highly incomplete frequency information. Information Theory, IEEE Transactionson, 52(2):489–509, 2006.
[2] E.J. Candes and T. Tao. Near-optimal signal recovery from random projections: Universalencoding strategies? Information Theory, IEEE Transactions on, 52(12):5406–5425, 2006.
[3] D.L. Donoho. Compressed sensing. Information Theory, IEEE Transactions on,52(4):1289–1306, 2006.
[4] A. Bilgin, M.W. Marcellin, and M.I. Altbach. Compression of electrocardiogram signals usingJPEG2000. Consumer Electronics, IEEE Transactions on, 49(4):833–840, 2003.
[5] D.L. Donoho and X. Huo. Uncertainty principles and ideal atomic decomposition. InformationTheory, IEEE Transactions on, 47(7):2845–2862, 2001.
[6] R. Coifman, F. Geshwind, and Y. Meyer. Noiselets. Applied and Computational HarmonicAnalysis, 10(1):27–44, 2001.
[7] J.F. Claerbout and F. Muir. Robust modeling with erratic data. Geophysics, 38:826, 1973.
[8] F. Santosa and W.W. Symes. Linear inversion of band-limited reflection seismograms. SIAMJournal on Scientific and Statistical Computing, 7:1307, 1986.
[9] J. Tropp and A.C. Gilbert. Signal recovery from partial information via orthogonal matchingpursuit, 2005.
[10] E. Cand`es and J. Romberg. Sparsity and incoherence in compressive sampling. Inverseproblems, 23:969, 2007.
[11] A.C. Gilbert, S. Muthukrishnan, and M. Strauss. Improved time bounds for near-optimalsparse Fourier representations. In Proceedings of SPIE, volume 5914, page 59141A. Citeseer,2005.
[12] M. Vetterli, P. Marziliano, and T. Blu. Sampling signals with finite rate of innovation. SignalProcessing, IEEE Transactions on, 50(6):1417–1428, 2002.
[13] D.L. Donoho and P.B. Stark. Uncertainty principles and signal recovery. SIAM Journal onApplied Mathematics, 49(3):906–931, 1989.
[14] P. Feng and Y. Bresler. Spectrum-blind minimum-rate sampling and reconstruction of multibandsignals. In icassp, pages 1688–1691. IEEE, 1996.
[15] E.J. Candes and T. Tao. Decoding by linear programming. Information Theory, IEEETransactions on, 51(12):4203–4215, 2005.
[16] E.J. Candes, J.K. Romberg, and T. Tao. Stable signal recovery from incomplete and inaccuratemeasurements. Communications on Pure and Applied Mathematics, 59(8):1207–1223,2006.
[17] EJ Cand`es. Lectures on compressive sampling and frontiers in signal processing. The Institutefor Mathematics and its Applications, pages 2006–2007.
[18] A. Cohen, W. Dahmen, and R. DeVore. Compressed sensing and best k-term approximation.American Mathematical Society, 22(1):211–231, 2009.
[19] E. Candes and T. Tao. The Dantzig selector: Statistical estimation when p is much largerthan n. The Annals of Statistics, 35(6):2313–2351, 2007.
[20] J. Haupt and R. Nowak. Signal reconstruction from noisy random projections. InformationTheory, IEEE Transactions on, 52(9):4036–4048, 2006.
[21] R. Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal StatisticalSociety. Series B (Methodological), 58(1):267–288, 1996.
[22] S.S. Chen, D.L. Donoho, and M.A. Saunders. Atomic decomposition by basis pursuit. SIAMjournal on scientific computing, 20(1):33–61, 1999.
[23] R. Baraniuk, M. Davenport, R. DeVore, and M. Wakin. A simple proof of the restrictedisometry property for random matrices. Constructive Approximation, 28(3):253–263, 2008.
[24] S. Mendelson, A. Pajor, and N. Tomczak-Jaegermann. Uniform uncertainty principle forBernoulli and subgaussian ensembles. Constructive Approximation, 28(3):277–289, 2008.
[25] M. Rudelson and R. Vershynin. On sparse reconstruction from Fourier and Gaussian measurements.Communications on Pure and Applied Mathematics, 61(8):1025–1045, 2008.
[26] R.E. Blahut. Algebraic codes for data transmission. Cambridge Univ Pr, 2003.
[27] D. Baron, M.B. Wakin, M.F. Duarte, S. Sarvotham, and R.G. Baraniuk. Distributed compressedsensing. preprint, pages 1–50, 2005.
[28] M. Lustig, D.L. Donoho, and J.M. Pauly. Rapid MR imaging with compressed sensing andrandomly under-sampled 3DFT trajectories. In Proc. 14th Ann. Meeting ISMRM. Citeseer.
[29] D. Takhar, V. Bansal, M. Wakin, M. Duarte, D. Baron, KF Kelly, and RG Baraniuk. Acompressed sensing camera: New theory and an implementation using digital micromirrors.Proc. Computational Imaging IV at SPIE Electronic Imaging, San Jose, 2006.
[30] E.J. Candes, M.B. Wakin, and S.P. Boyd. Enhancing sparsity by reweighted l1 minimization.Journal of Fourier Analysis and Applications, 14(5):877–905, 2008.


2011-03-16

Fast Multidimensional Scaling using Vector Extrapolation

本文Google Docs地址

Guy Rosman,Alexander M. Bronstein,Michael M. Bronstein,Avram Sidi,Ron Kimmel
Fast Multidimensional Scaling using Vector Extrapolation
Technion - Computer Science Department - Technical Report CIS-2008-01 --2008

摘要

MDS是一类低维表示方法,对象是给定距离矩阵的点集。在很多MDS应用中,算法的速

度和效率是关键问题,向量外推方法可以用来提高固定点迭代算法的收敛速度。本文

将提出用向量外推加速MDS数值解速度的方法。

2 多维标度

2.1 最小二乘MDS

在n维欧式空间中,坐标表达:$latex {X=\left( {x_{ij} } \right)}&fg=000000$

,是N$latex {\times }&fg=000000$m度量值;$latex {d_{ij} \left( X \right)}

&fg=000000$表示i、j之间的距离;$latex {\delta _{ij} }&fg=000000$表示输入距

离,本文将使用一个新的表示形式。

通常MDS问题可以用下述优化模型表示:

$latex \displaystyle


\mathop {\min }\limits_X \sum\limits_{i<j} {w_{ij} f_{ERR} \left( {d_

{ij} \left( X \right),\delta _{ij} } \right)} &fg=000000$

$latex {w_{ij} }&fg=000000$是表征每一对距离重要性的权值;$latex {f_{ERR}

}&fg=000000$是计算输入距离和近似距离的误差代价函数,选择距离误差函数的二次

项形式,即「stress 」函数:

$latex \displaystyle


f_{STRESS} \left( {d,\delta } \right)=\left( {d-\delta } \right)^2

&fg=000000$



定义距离平方之差的形式:

$latex \displaystyle f_{SSTRESS}


\left( {d,\delta } \right)=\left( {d^2-\delta ^2} \right)^2 &fg=000000$

得到一个函数,通常称之为「sstress」函数。

2.2 SMACOF 算法

为了最小化stress函数,$latex {s\left( {\rm

{\bf X}} \right)=\sum\limits_{i<j}^ {w_{ij} \left( {d_{ij} \left( {\rm

{\bf X}} \right)-\delta _{ij} } \right)^2} }&fg=000000$,我们需要求出$latex

{s\left( {\rm {\bf X}} \right)}&fg=000000$关于$latex {{\rm {\bf X}}}

&fg=000000$梯度:

$latex \displaystyle \nabla s\left( {\rm


{\bf X}} \right)=2{\rm {\bf VX}}-2{\rm {\bf B}}\left( {\rm {\bf X}}

\right){\rm {\bf X}} &fg=000000$



这里$latex {{\rm {\bf V}}}&fg=000000$和$latex {{\rm {\bf B}}}&fg=000000$通

过下式给出:

$latex \displaystyle \left(




{\rm {\bf V}} \right)_{ij} =\left\{ {\begin{array}{c} -w_{ij} \mbox{ } if



i\ne j \\ \sum\limits_{k\ne i}^ {w_{ik} } if i=j \\ \end{array}} \right. \



\ \ \ \ (1)&fg=000000$




$latex \displaystyle \left( {\rm {\bf




B}} \right)_{ij} =\left\{ {\begin{array}{c} \mbox{-w}_{ij} \delta _{ij} d_



{ij}^{-1} \left( {\rm {\bf X}} \right) \mbox{if }i\ne j \mbox{and} d_{ij}



\left( {\rm {\bf X}} \right)\ne 0 \\ 0 \mbox{if }i\ne j \mbox{and} d_{ij}



\left( {\rm {\bf X}} \right)=0 \\ -\sum\nolimits_{k\ne i} {b_{ik} } \mbox



{if }i=j \\ \end{array}} \right. \ \ \ \ \ (2)&fg=000000$



使用一阶优化,则通过下式得到:

$latex {2{\rm {\bf VX}}=2{\rm {\bf B}}\left( {\rm {\bf X}} \right){\rm {\bf

X}}}&fg=000000$,或者:

$latex \displaystyle




{\rm {\bf X}}={\rm {\bf V}}^\dag {\rm {\bf B}}\left( {\rm {\bf X}}



\right){\rm {\bf X}} \ \ \ \ \ (3)&fg=000000$



这里$latex {\dag }&fg=000000$表示矩阵伪逆。

由(3)可以得到迭代形式(方法见文献[25]):

$latex \displaystyle {\rm {\bf X}}^{\left( {k+1} \right)}={\rm {\bf V}}^

\dag {\rm {\bf B}}\left( {{\rm {\bf X}}^{\left( k \right)}} \right){\rm

{\bf X}}^{\left( k \right)} \ \ \ \ \ (4)&fg=000000$

(3)式可视为一个固定点,将(3)迭代使之收敛到stress 代价函数的局部极

小值。以上处理思路称为「SMACOF」------standing from Scaling by Majorizing a

Complicated Function。

stress函数有一个重要特性就是可以保证一个stress值的单调下降序列。这是其他优

化问题少见的。另一方面,SMACOF算法的收敛速度是比较慢的。

2.3 经典标度

另一类使用广泛的代价函数是「strain

(定义见下文),其引出了一类代数MDS算法,称为「经典标度」(classical

scaling)。

令$latex {\Delta ^{\left( 2 \right)}}&fg=000000$表示输入距离矩阵的平方;

$latex {{\rm {\bf D}}^{\left( \mbox{2} \right)}\left( {\rm {\bf X}}

\right)}&fg=000000$表示目标欧式距离的平方:

$latex \displaystyle {\rm {\bf D}}^\mbox{2}\left( {\rm {\bf X}}\right)\mbox{=}{\rm {\bf c1}}^T+{\rm {\bf 1c}}^T-2{\rm {\bf XX}}^T \ \ \ \ \ (5)&fg=000000$

这里$latex {c_i =\left\langle {x_i ,x_i } \right\rangle }&fg=000000$

令$latex {{\rm {\bf J}}\mbox{=}{\rm {\bf I}}\mbox{-}\frac{1}{N}{\rm {\bf

11}}^T}&fg=000000$,表示中心矩阵。

假设给定距离均为欧式距离,得到这些点的内积矩阵(Gram 矩阵):

$latex \displaystyle {\rm {\bf B}}_\Delta =\frac{1}{2}{\rm


{\bf J}}\Delta ^{\left( 2 \right)}{\rm {\bf J}} &fg=000000$

对于欧式距离的情形,$latex {\left( {{\rm {\bf B}}_\Delta } \right)_{ij} =

\left\langle {x_i ,x_j } \right\rangle }&fg=000000$

但实际上,输入矩阵往往不是欧式的,二者的差异度可以通过一个测量值近似表示,

称之为「strain」:

$latex \displaystyle \begin


{array}{l} \left\| {-\frac{1}{2}{\rm {\bf J}}\left( {{\rm {\bf D}}^{\left(

2 \right)}\left( {\rm {\bf X}} \right)-\Delta ^{\left( 2 \right)}} \right)

{\rm {\bf J}}} \right\|_F^2 \\ {\rm {\bf =}}\left\| {{\rm {\bf XX}}^T+

\frac{1}{2}{\rm {\bf J}}\left( {\Delta ^{\left( 2 \right)}} \right){\rm

{\bf J}}} \right\|_F^2 \\ =\left\| {{\rm {\bf XX}}^T-{\rm {\bf B}}_\Delta }

\right\|_F^2 \\ \end{array} &fg=000000$

通过对$latex {{\rm {\bf B}}_\Delta }&fg=000000$进行特征值分解,可以找到

strain函数的全局最优解。

经典标度方法的缺陷在于它的自适应性不强,而且计算开销也较大。因此在后文的论

述中,将只关注stress函数以及SMACOF算法。然而,全局优化可保证经典标度方法在

SMACOF的初始阶段具有良好性能。

3 向量外推方法

考虑使用历次迭代方法来加速SMACOF算法的收敛速度,如使用向量外推方法来预测收

敛极限。

本文考察了两个向量外推方法:

1、 「MPE」------最小多项式外推(minimal polynomial extrapolation,文献[11]

);

2、「RRE」------减秩外推(reduced rank extrapolation,文献[31,19])。

两个方法在加速非线性/线性/大型稀疏系统方程中固定点迭代的向量序列收敛速度方

面都具有高效的性能。

二者均需要考察一个经过线性固定点迭代过程得出的序列($latex {{\rm {\bf x}}_0

,{\rm {\bf x}}_1 ,...}&fg=000000$):

$latex




\displaystyle x_{n+1} ={\rm {\bf A}}x_n +b, n=0,1,... \ \ \ \ \ (6)



&fg=000000$



此处$latex {{\rm {\bf A}}}&fg=000000$是给定$latex {N\times N}

&fg=000000$矩阵,$latex {{\rm {\bf b}}}&fg=000000$是给定N维向量,$latex

{x_0 }&fg=000000$是用户选取的初始向量。这个序列有一个极限s,它是下述方程的

唯一解:

$latex \displaystyle {\rm {\bf




x}}={\rm {\bf Ax}}+{\rm {\bf b}} \ \ \ \ \ (7)&fg=000000$



假设$latex {\rho \left( {\rm {\bf A}} \right)}&fg=000000$是谱半

径,$latex {\rho \left( {\rm {\bf A}} \right)<{\rm x}}&fg=000000$,另一

方面(7)式也可以写为$latex {\left( {{\rm {\bf I}}-{\rm {\bf A}}} \right)

x={\rm {\bf b}}}&fg=000000$,由于1不是$latex {{\rm {\bf A}}}&fg=000000$的奇

异值,所以矩阵$latex {{\rm {\bf I}}-{\rm {\bf A}}}&fg=000000$是非奇异的。

给定(6)式的序列,令

$latex \displaystyle


{\rm {\bf u}}_{\rm {\bf n}} {\rm {\bf =\Delta x}}_{\rm {\bf n}} {\rm {\bf

=x}}_{{\rm {\bf n+1}}} {\rm {\bf -x}}_{\rm {\bf n}} {\rm {\bf , n=0,1,...}}

&fg=000000$

再定义误差向量:

$latex \displaystyle ?_n




=x_n -s\mbox{, }n=0,1,... \ \ \ \ \ (8)&fg=000000$



考虑$latex {{\rm {\bf s}}={\rm {\bf As}}+{\rm {\bf b}}}&fg=000000$,我

们可以通过n步迭代从初始误差得到当前误差:

$latex \displaystyle ?_n =\left( {{\rm {\bf Ax}}_{n-1} +b} \right)-\left(

{{\rm {\bf As}}+{\rm {\bf b}}} \right)={\rm {\bf A}}\left( {{\rm {\bf x}}_

{n-1} -{\rm {\bf s}}} \right)={\rm {\bf A}}?_{n-1} \ \ \ \ \ (9)

&fg=000000$

由此可以得到:

$latex \displaystyle




?_n ={\rm {\bf A}}^n?_0 \mbox{, }n=0,1,... \ \ \ \ \ (10)&fg=000000$



通过对k+1连续$latex {x_i }&fg=000000$(k +1 consecutive)采用「加权平

均」来近似s:

$latex \displaystyle




{\rm {\bf s}}_{n,k} =\sum\limits_{i=0}^k {\gamma _i {\rm {\bf x}}_{n+i} }



\mbox{; }\sum\limits_{i=0}^k {\gamma _i } =1 \ \ \ \ \ (11)&fg=000000$



将(8)带入(11),并考虑$latex {\sum\nolimits_{i=0}

^k {\gamma _i } =1}&fg=000000$,可以得到:

$latex \displaystyle {\rm {\bf s}}_{n,k} =\sum\limits_{i=0}^k {\gamma _i

\left( {{\rm {\bf s}}+?_{n+i} } \right)} ={\rm {\bf s}}+\sum\limits_{i=0}^k

{\gamma _i } ?_{n+i} \ \ \ \ \ (12)&fg=000000$

再带入(10),得到

$latex \displaystyle {\rm {\bf s}}_{n,k} ={\rm {\bf s}}+\sum\limits_

{i=0}^k {\gamma _i {\rm {\bf A}}^{n+i}} ?_0 \ \ \ \ \ (13)&fg=000000$

为了使$latex {?_{n+i\mbox{, }} i=0,1,...}&fg=000000$的加权和$latex

{\sum\nolimits_{i=0}^k {\gamma _i {\rm {\bf A}}^{n+i}?_0 } }&fg=000000$尽可

能小,我们必须选择合适的$latex {\gamma _i }&fg=000000$。

现在,给定一个N$latex {\times }&fg=000000$N矩阵B,以及任意一个N维向

u,则存在具有最小阶数(最多N)的唯一莫尼多项式(monic polynomial)

$latex {P\left( z \right)}&fg=000000$,是u的零化多项式:$latex {P

\left( {\rm {\bf B}} \right){\rm {\bf u}}=0}&fg=000000$,这个多项式称为

B关于u的「极小多项式」(minimal polynomial)。$latex {P\left(

z \right)}&fg=000000$的零解部分或者全部是B的特征值。

因此,若A关于$latex {?_n }&fg=000000$的极小多项式是:

$latex \displaystyle P\left( z \right)=\sum\limits_{i=0}^k


{c_i z^i} ;\mbox{ }c_k =1 &fg=000000$

即:$latex {P\left( {\rm {\bf A}} \right)?_n =0}&fg=000000$,联立(10)式,得:

$latex




\displaystyle \sum\limits_{i=0}^k {c_i {\rm {\bf A}}^i?_n } =\sum



\limits_{i=0}^k {c_i ?_{n+i} } =0 \ \ \ \ \ (14)&fg=000000$



(14)是N个线性方程的集,k个未知参数$latex {c_0 ,c_1

,...,}&fg=000000$同时$latex {c_k =1}&fg=000000$。这N个方程的解是唯一的,因

为$latex {P\left( z \right)}&fg=000000$的解是唯一的。我们可以由x$latex {_

{i}}&fg=000000$的信息单独得到c$latex {_{i}}&fg=000000$,联立(14)和(9),有:

$latex


\displaystyle 0=\sum\limits_{i=0}^k {c_i {\rm {\bf A}}?_{n+i} } =\sum

\limits_{i=0}^k {c_i ?_{n+i+1} } &fg=000000$

减去(14),由此得到线性系统:

$latex \displaystyle \sum\limits_{i=0}^k {c_i {\rm {\bf u}}




_{n+i} } ={\rm {\bf 0}} \ \ \ \ \ (15)&fg=000000$



一旦$latex {c_0 ,c_1 ,...,c_{k-1} }&fg=000000$确定,我们令$latex {c_k

}&fg=000000$=1、$latex {\gamma _i ={c_i } \mathord{\left/ {\vphantom {{c_i

} {\sum\nolimits_{j=0}^k {c_j } }}} \right. \kern-\nulldelimiterspace}

{\sum\nolimits_{j=0}^k {c_j } },\mbox{ }i=0,1,...,k.}&fg=000000$

综上所述,若k是A关于$latex {?_n }&fg=000000$的极小多项式的阶数,则存在序列

满足$latex {\sum\nolimits_{i=0}^k {\gamma _i =1} }&fg=000000$,则有:$latex

{\sum\nolimits_{i=0}^k {\gamma _i {\rm {\bf x}}_{n+i} ={\rm {\bf s}}} }

&fg=000000$(Why not $latex {{\rm {\bf s}}_{n,k} }&fg=000000$?)

值得注意的是,不论是否$latex {\rho \left( {\rm {\bf A}} \right)<1}

&fg=000000$,s都是$latex {\left( {{\rm {\bf I}}-{\rm {\bf A}}}

\right){\rm {\bf x}}={\rm {\bf b}}}&fg=000000$的解,因此,无论$latex

{\mathop {\lim }\limits_{n\rightarrow \infty } {\rm {\bf x}}_n }&fg=000000$

是否存在,$latex {{\rm {\bf s}}=\sum\nolimits_{i=0}^k {\gamma _i {\rm {\bf

x}}_{n+i} } }&fg=000000$都成立。

为简化表述,使用如下标记:

$latex




\displaystyle {\rm {\bf U}}_s^{\left( j \right)} =\left[ {{\rm {\bf u}}_j



\vert {\rm {\bf u}}_{j+1} \vert ...\vert {\rm {\bf u}}_{j+s} } \right] \ \



\ \ \ (16)&fg=000000$



因此,$latex {{\rm {\bf U}}_s^{\left( j \right)} }&fg=000000$是一个N

$latex {\times }&fg=000000$(j+1)矩阵,则(14)可以表述

为:

$latex \displaystyle {\rm {\bf U}}




_k^{\left( n \right)} {\rm {\bf c}}={\rm {\bf 0}}\mbox{; }{\rm {\bf c}}=



\left[ {c_0 ,c_1 ,...c_k } \right]^T \ \ \ \ \ (17)&fg=000000$



当然,用$latex {\sum\nolimits_{i=0}^k {c_i } }&fg=000000$分解(17),

可以得到:

$latex \displaystyle {\rm {\bf




U}}_k^{\left( n \right)} \gamma ={\rm {\bf 0}}\mbox{; }\gamma =\left[



{\gamma _0 ,\gamma _1 ,...\gamma _k } \right]^T \ \ \ \ \ (18)



&fg=000000$





3.1 MPE派生

前文所述,极小多项式的阶数可以达到N,若N很大

,这就可能会导致比较大的存储开销;此外,我们还没有方法可以准确得知这个阶数

,鉴于此,可以采取一种解决途径:先任意找一个远远小于阶数的正指数k,带入(15),则(15)则不再是一致的,因此也不再

有对于$latex {c_0 ,c_1 ,...,c_{k-1} }&fg=000000$(其中$latex {c_\mbox{k}

=1}&fg=000000$)的唯一解。通常对于这样的问题,我们可以采取最小二乘方法求解

,沿着这个思路,计算(15)中描述的$latex {\gamma _0 ,

\gamma _1 ,...\gamma _k }&fg=000000$,然后计算向量$latex {{\rm {\bf s}}_

{n,k} =\sum\limits_{i=0}^k {\gamma _i {\rm {\bf x}}_{n+i} } }&fg=000000$,

计算结果即为s的近似。以上方法称为「极小多项式外推算法」(MPE------

minimal polynomial extrapolation),其算法流程见表-1:

\centerline{\includegraphics[width=5.83in,height=1.41in]{mytex31.eps}}

\caption{MPE算法流程}

3.2 RRE派生

同3.1,在(18)中带入k,则不再对$latex {\gamma _0 ,

\gamma _1 ,...\gamma _k }&fg=000000$具有唯一解,同样应用最小二乘算法(约束

为$latex {\sum\nolimits_{i=0}^k {\gamma _i =1} }&fg=000000$),然后计算

s的近似$latex {{\rm {\bf s}}_{n,k} =\sum\limits_{i=0}^k {\gamma _i

{\rm {\bf x}}_{n+i} } }&fg=000000$。这个方法称为「减秩外推」(RRE------

reduced rank extrapolation),算法流程见表-2:

\centerline{\includegraphics[width=5.83in,height=1.22in]{mytex32.eps}}

\caption{RRE算法}

3.3 对非线性方程的处理

前文所述,在SMACOF算法中存在非线性方程,现在用向量外推方法处理这个问题。假

设非线性方程记作:

$latex \displaystyle




{\rm {\bf x}}={\rm {\bf F}}\left( {\rm {\bf x}} \right) \ \ \ \ \ (19)



&fg=000000$



这里$latex {{\rm {\bf F}}\left( {\rm {\bf x}} \right)}&fg=000000$是N维

向量值函数,x是N维未知向量,通过下式得到近似值x$latex {_{n}}

&fg=000000$:

$latex \displaystyle {\rm




{\bf x}}_{n+1} ={\rm {\bf F}}\left( {{\rm {\bf x}}_n } \right)\mbox{, }



n=0,1,... \ \ \ \ \ (20)&fg=000000$



并且假设这个序列收敛于解s,$latex {{\rm {\bf F}}}&fg=000000$是

SMACOF迭代的右边函数,由于x接近s,$latex {{\rm {\bf F}}\left(

{\rm {\bf x}} \right)}&fg=000000$可通过泰勒级数展开:

$latex \displaystyle {\rm {\bf F}}\left( {\rm {\bf x}} \right)={\rm {\bf

F}}\left( {\rm {\bf s}} \right)+{\rm {\bf {F}'}}\left( {\rm {\bf s}}

\right)\left( {{\rm {\bf x}}-{\rm {\bf s}}} \right)+o\left( {\left\| {{\rm

{\bf x}}-{\rm {\bf s}}} \right\|^2} \right)\mbox{ as }{\rm {\bf

x}}\rightarrow {\rm {\bf s}} &fg=000000$

此处$latex {{\rm {\bf {F}'}}\left( {\rm {\bf x}} \right)}&fg=000000$是

$latex {{\rm {\bf F}}\left( {\rm {\bf x}} \right)}&fg=000000$的雅各比矩阵,

又因为$latex {{\rm {\bf F}}\left( {\rm {\bf s}} \right)={\rm {\bf s}}}

&fg=000000$,则上式可写为:

$latex \displaystyle {\rm {\bf


F}}\left( {\rm {\bf x}} \right)={\rm {\bf s}}+{\rm {\bf {F}'}}\left( {\rm

{\bf s}} \right)\left( {{\rm {\bf x}}-{\rm {\bf s}}} \right)+o\left(

{\left\| {{\rm {\bf x}}-{\rm {\bf s}}} \right\|^2} \right)\mbox{ as }{\rm

{\bf x}}\rightarrow {\rm {\bf s}} &fg=000000$

假设序列$latex {{\rm {\bf x}}_0 ,{\rm {\bf x}}_1 ,...,}&fg=000000$收敛于

s,则当n足够大时,$latex {{\rm {\bf x}}_n }&fg=000000$逼近s

因此有:

$latex \displaystyle {\rm {\bf x}}_{n+1} ={\rm


{\bf s}}+{\rm {\bf {F}'}}\left( {\rm {\bf s}} \right)\left( {{\rm {\bf x}}

_n -{\rm {\bf s}}} \right)+o\left( {\left\| {{\rm {\bf x}}_n -{\rm {\bf

s}}} \right\|^2} \right)\mbox{ as n}\rightarrow \infty &fg=000000$

即:$latex {{\rm {\bf x}}_{n+1} -{\rm {\bf s}}={\rm {\bf {F}'}}\left( {\rm

{\bf s}} \right)\left( {{\rm {\bf x}}_n -{\rm {\bf s}}} \right)+o\left(

{\left\| {{\rm {\bf x}}_n -{\rm {\bf s}}} \right\|^2} \right)\mbox{ as

n}\rightarrow \infty }&fg=000000$。

对于所有的大n,向量$latex {{\rm {\bf x}}_n }&fg=000000$可视为从形如$latex

{\left( {{\rm {\bf I}}-{\rm {\bf A}}} \right){\rm {\bf x}}={\rm {\bf b}}}

&fg=000000$的线性系统产生:

$latex




\displaystyle {\rm {\bf x}}_{n+1} ={\rm {\bf Ax}}_n +{\rm {\bf b}}\mbox{,



}n=0,1,..., \ \ \ \ \ (21)&fg=000000$



此处$latex {{\rm {\bf A}}={\rm {\bf {F}'}}\left( {\rm {\bf s}}

\right)}&fg=000000$,$latex {{\rm {\bf b}}=\left[ {{\rm {\bf I}}-{\rm {\bf

{F}'}}\left( {\rm {\bf s}} \right)} \right]{\rm {\bf s}}}&fg=000000$。

MPE和RRE已经应用在大型稀疏数值系统中,如计算流体动力学、半导体研究和X-射线

断层扫描系统。

3.4 MPE和RRE的有效应用

MPE和RRE的关键问题是关于最小二乘问题的精确解

和尽可能减少计算时间和存储空间开销。

关于最小二乘问题的解决,可以采用对基$latex {{\rm {\bf U}}_k^{\left( n

\right)} }&fg=000000$进行QR分解:

$latex \displaystyle


{\rm {\bf U}}_k^{\left( n \right)} ={\rm {\bf Q}}_k {\rm {\bf R}}_k

&fg=000000$

此处$latex {{\rm {\bf Q}}_k }&fg=000000$是N$latex {\times }&fg=000000$(k

+1)阶酉矩阵,可记为

$latex \displaystyle




{\rm {\bf Q}}_k =\left[ {{\rm {\bf q}}_0 \vert {\rm {\bf q}}_1 \vert ...



\vert {\rm {\bf q}}_k } \right] \ \ \ \ \ (22)&fg=000000$



其中,$latex {{\rm {\bf q}}_i^\ast {\rm {\bf q}}_j =\delta _{ij} }

&fg=000000$;

$latex {{\rm {\bf R}}_k }&fg=000000$是(k+1)$latex {\times }

&fg=000000$(k+1)阶上三角矩阵,对角线元素均为正数:

$latex \displaystyle {\rm {\bf R}}_k =\left[ {{\begin




{array}{*{20}c} {r_{00} } \hfill } \hfill & \cdots \hfill } \hfill \\ \hfill } \hfill & \cdots \hfill }



\hfill \\ \hfill & \hfill & \ddots \hfill & \vdots \hfill \\ \hfill &



\hfill & \hfill } \hfill \\ \end{array} }} \right]\mbox{; }r_{ii}



\mbox{>0, }i=0,1,...,k \ \ \ \ \ (23)&fg=000000$



QR分解可采用改良Gram-Schmidt正交方法(MGS),表-3描述了对$latex {{\rm

{\bf U}}_k^{\left( n \right)} }&fg=000000$应用MGS的流程:

\centerline{\includegraphics[width=5.93in,height=1.48in]{mytex33.eps}}

\caption{MGS算法}

此处$latex {\left\| {\rm {\bf x}} \right\|\mbox{-}\sqrt {{\rm {\bf x}}^\ast

{\rm {\bf x}}} }&fg=000000$,$latex {{\rm {\bf u}}_i^{\left( {j+1} \right)}

}&fg=000000$重写为$latex {{\rm {\bf u}}_i^{\left( j \right)} }&fg=000000$,

如此则$latex {{\rm {\bf u}}_{n+i} ,{\rm {\bf u}}_i^{\left( j \right)} ,{\rm

{\bf q}}_i }&fg=000000$占有相同的存储空间。

很重要的一点是,当构建矩阵$latex {{\rm {\bf U}}_k^{\left( n \right)} }

&fg=000000$时,我们将$latex {{\rm {\bf x}}_{n+i} }&fg=000000$重写为$latex

{{\rm {\bf u}}_{n+i} ={\rm {\bf x}}_{n+i} }&fg=000000$,而只保存$latex

{{\rm {\bf x}}_n }&fg=000000$;下一步,当计算矩阵$latex {{\rm {\bf Q}}_k }

&fg=000000$时,我们把$latex {{\rm {\bf q}}_i \mbox{,i=0,1,...,k}}

&fg=000000$重写为$latex {{\rm {\bf u}}_{n+i} }&fg=000000$。这就意味着,在计

算$latex {{\rm {\bf Q}}_k }&fg=000000$和$latex {{\rm {\bf R}}_x }

&fg=000000$的每一阶段,我们都只保留k+2个向量,$latex {{\rm {\bf x}}_{n+1}

,...,{\rm {\bf x}}_{n+k+1} }&fg=000000$无需保存。

表-4展示了采用QR分解的MPE和RRE算法,它们具有一致的框架:

\centerline{\includegraphics[width=6.00in,height=5.63in]{mytex34.eps}}

\caption{采用QR分解的MPR/RRE算法}

3.5 误差估计

1、 对于线性序列:当(6)中的迭代向量$latex {{\rm {\bf

x}}_i }&fg=000000$是线性时,有:

$latex \displaystyle {\rm


{\bf r}}\left( {\rm {\bf x}} \right)={\rm {\bf b}}-\left( {{\rm {\bf I}}-

{\rm {\bf A}}} \right){\rm {\bf x}}=\left( {{\rm {\bf Ax}}+{\rm {\bf b}}}

\right)-{\rm {\bf x}} &fg=000000$

即:$latex {{\rm {\bf r}}\left( {{\rm {\bf x}}_n } \right)={\rm {\bf x}}_

{n+1} -{\rm {\bf x}}_n ={\rm {\bf u}}_n }&fg=000000$

在应用MPE和RRE算法时,考虑$latex {\sum\nolimits_{i=0}^k {\gamma _i =1} }

&fg=000000$,可得

$latex \displaystyle




{\rm {\bf r}}\left( {{\rm {\bf s}}_{n,k} } \right)=\sum\limits_{i=0}^k



{\gamma _i {\rm {\bf u}}_{n+i} } ={\rm {\bf U}}_k^{\left( n \right)} \gamma



\ \ \ \ \ (24)&fg=000000$



我们考察其$latex {l_2 }&fg=000000$范数,即:

$latex


\displaystyle \left\| {{\rm {\bf r}}\left( {{\rm {\bf s}}_{n,k} } \right)}

\right\|=\left\| {{\rm {\bf U}}_k^{\left( n \right)} \gamma } \right\|

&fg=000000$

2、 对于非线性序列:当(20)中的$latex {{\rm {\bf x}}_i

}&fg=000000$是线性时,有:

$latex \displaystyle {\rm {\bf


r}}\left( {{\rm {\bf s}}_{n,k} } \right)={\rm {\bf F}}\left( {{\rm {\bf

s}}_{n,k} } \right)-{\rm {\bf s}}_{n,k} \approx {\rm {\bf U}}_k^{\left( n

\right)} \gamma &fg=000000$

因此:

$latex \displaystyle \left\| {{\rm {\bf r}}\left(


{{\rm {\bf s}}_{n,k} } \right)} \right\|\approx \left\| {{\rm {\bf U}}_k^

{\left( n \right)} \gamma } \right\| &fg=000000$

不论$latex {{\rm {\bf x}}_i }&fg=000000$是否是线性的,$latex {\left\|

{{\rm {\bf U}}_k^{\left( n \right)} \gamma } \right\|}&fg=000000$都可以无需

计算$latex {{\rm {\bf s}}_{n,k} }&fg=000000$而得出:

$latex


\displaystyle \left\| {{\rm {\bf U}}_k^{\left( n \right)} \gamma } \right

\|=\left\{ {{\begin{array}{*{20}c} {r_{kk} \left| {\gamma _k } \right|

\mbox{ for MPE}} \hfill \\ {\sqrt \lambda \mbox{ for RRE}} \hfill \\ \end

{array} }} \right. &fg=000000$

此处$latex {r_{kk} }&fg=000000$是矩阵$latex {{\rm {\bf R}}_k }&fg=000000$

对角线上的最后一个元素,$latex {\lambda }&fg=000000$是上一节算法的第三步中

得到的参数(文献[44])。

3.6 MPE和RRE的误差分析

关于MPE和RRE的线性误差分析在文献[42、46、45、48、49]中已有论述。本文将重点

介绍非线性的情况下,向量外推方法可以达到怎样的性能。

定理1 假设向量$latex {{\rm {\bf x}}_n }&fg=000000$满足:

$latex \displaystyle {\rm {\bf x}}_n ={\rm {\bf s}}+\sum


\limits_{i=1}^p {{\rm {\bf v}}_i \lambda _i^n } &fg=000000$

此处,向量$latex {{\rm {\bf v}}_i }&fg=000000$是线性独立的,非零标量$latex

{\lambda _i \ne 1}&fg=000000$,取值不同,且满足:

$latex


\displaystyle \left| {\lambda _1 } \right|\ge \left| {\lambda _2 }

\right|\ge ... &fg=000000$

若有$latex {\left| {\lambda _k } \right|\ge \left| {\lambda _{k+1} }

\right|}&fg=000000$,则对于MPE和RRE,有:

$latex


\displaystyle {\rm {\bf s}}_{n,k} -{\rm {\bf s}}=o\left( {\lambda _{k+1}^n

} \right)\mbox{ as }n\rightarrow \infty &fg=000000$

以及:

$latex \displaystyle \mathop {\lim }\limits_{n


\rightarrow \infty } \sum\limits_{i=0}^k {\gamma _i^{\left( {n,k} \right)}

z^i} =\prod\limits_{i=1}^k {\frac{\lambda -\lambda _i }{1-\lambda _i }}

&fg=000000$

此处,$latex {\gamma _i^{\left( {n,k} \right)} }&fg=000000$代表$latex

{\gamma _i }&fg=000000$。

当$latex {{\rm {\bf x}}_n }&fg=000000$是(6)中迭代产生时

,$latex {\lambda _i }&fg=000000$部分或者全部是迭代矩阵$latex {{\rm {\bf

A}}}&fg=000000$($latex {{\rm {\bf A}}}&fg=000000$是对角化的)的非零特征值

,$latex {{\rm {\bf v}}_i }&fg=000000$是对应的特征向量。

定理2 假设向量$latex {{\rm {\bf x}}_n }&fg=000000$由$latex {{\rm

{\bf x}}_{n+1} ={\rm {\bf Ax}}_n +{\rm {\bf b}}}&fg=000000$迭代产生,矩阵

$latex {{\rm {\bf I}}-{\rm {\bf A}}}&fg=000000$非奇异,$latex {\left( {{\rm

{\bf I}}-{\rm {\bf A}}} \right){\rm {\bf x}}={\rm {\bf b}}}&fg=000000$的解

为$latex {{\rm {\bf s}}}&fg=000000$,$latex {{\rm {\bf A}}}&fg=000000$的非

零特征值为:

$latex \displaystyle \left| {\lambda _1 }


\right|\ge \left| {\lambda _2 } \right|\ge ... &fg=000000$

不论是否$latex {\left| {\lambda _k } \right|\ge \left| {\lambda _{k+1} }

\right|}&fg=000000$,对于

(i)RRE无条件满足;

(ii)MPE在提供$latex {{\rm {\bf I}}-{\rm {\bf A}}}&fg=000000$的特征值(这

些特征值均位于一条复平面上通过原点的直线的同一侧,例如当$latex {{\rm {\bf

A}}+{\rm {\bf A}}^\ast }&fg=000000$是正定)时,

下式均成立:

$latex \displaystyle {\rm {\bf s}}_{n,k} -


{\rm {\bf s}}=o\left( {\lambda _{k+1}^n } \right)\mbox{ as }n\rightarrow

\infty &fg=000000$

综合定理1、2,注意到我们并未假设$latex {\lim _{n\rightarrow \infty } {\rm

{\bf x}}_n }&fg=000000$存在,对于$latex {\lim _{n\rightarrow \infty } {\rm

{\bf x}}_n }&fg=000000$不存在的情况,$latex {\left| {\lambda _1 } \right|

\ge 1}&fg=000000$的条件是必须的。

事实上,若用$latex {{\rm {\bf x}}_i }&fg=000000$逼近s,可得误差

$latex \displaystyle ?_n ={\rm {\bf x}}_n -{\rm {\bf s}}=o


\left( {\lambda _1^n } \right)\mbox{ as }n\rightarrow \infty

&fg=000000$

因此,无论$latex {\lim _{n\rightarrow \infty } {\rm {\bf x}}_n }

&fg=000000$存在与否,只要给定$latex {\left| {\lambda _{k+1} } \right|

<1}&fg=000000$,我们均可得到$latex {\lim _{n\rightarrow \infty } {\rm

{\bf s}}_{n,k} ={\rm {\bf s}}}&fg=000000$。此外,当给定$latex {\left|

{\lambda _{k+1} } \right|<\left| {\lambda _1 } \right|}&fg=000000$时,序

列$latex {{\rm {\bf x}}_n }&fg=000000$收敛于s、$latex {{\rm {\bf

s}}_{n,k} }&fg=000000$收敛于s的速度更快。即,MPE和RRE方法加速了序列

$latex {\left\{ {{\rm {\bf x}}_n } \right\}}&fg=000000$的收敛速度。

3.7 循环MPE/RRE

定理1、2需要固定k,然后令n趋于无限大,很显然

实际中无法满足这个条件;此外增加k也可以使得MPE和RRE的收敛速度加快,然而我们

同样无法无限增大k。

循环(或者再启动)方法可以用来解决上述问题,在循环模式中,n和k是固定的,

表-5的算法流程中1~3步称为一个「循环」,第i次循环的$latex {{\rm {\bf s}}_

{n,k} }&fg=000000$记作$latex {{\rm {\bf s}}_{n,k}^{\left( i \right)} }

&fg=000000$。

\centerline{\includegraphics[width=5.91in,height=1.30in]{mytex35.eps}}

\caption{循环模式下的向量外推算法}

循环模式带来的两个好处(文献[48,49]):

1、产生一个紧致的误差上界;

2、防止了「广义最小残差停滞」(GMRES stagnates)

循环MPE/RRE在非线性系统下的性能分析(文献[50,51])带来的启发意义是,当第i次

循环中k接近k$latex {_{i}}&fg=000000$时,矩阵$latex {{\rm {\bf {F}'}}\left(

{\rm {\bf s}} \right)}&fg=000000$关于$latex {?_\mbox{0} \mbox{=}{\rm {\bf

x}}_\mbox{0} -{\rm {\bf s}}}&fg=000000$的极小多项式的阶数------序列$latex

{\left\{ {{\rm {\bf s}}_{n,k}^{\left( i \right)} } \right\}_{i=0}^\infty }

&fg=000000$将二阶收敛于$latex {{\rm {\bf s}}}&fg=000000$。

但由于k$latex {_{i}}&fg=000000$可以等于N,且无法知道它的准确大小,对于大规

模问题的应用,其存储开销可能会很大;另一方面,尝试从循环MPE/RRE实现二阶收敛

可能是不现实的。

3.8 与Krylov子空间方法结合

对于线性序列应用,MPE和RRE方法与Krylov子空间

方法------如Arnoldi(文献[1])、GMRES(文献[38])------关系很大。文献[43]给

出了如下定理:

定理3 考察线性系统$latex {\left( {{\rm {\bf I}}-{\rm {\bf A}}}

\right){\rm {\bf x}}={\rm {\bf b}}}&fg=000000$,其中$latex {{\rm {\bf x}}_0

}&fg=000000$是初始向量,令向量序列$latex {\left\{ {{\rm {\bf x}}_n }

\right\}}&fg=000000$通过$latex {{\rm {\bf x}}_{n+1} ={\rm {\bf Ax}}_n +{\rm

{\bf b}}}&fg=000000$生成。分别应用MPE、RRE,从该序列中生成$latex {{\rm {\bf

s}}_{0,k}^{MPE} }&fg=000000$和$latex {{\rm {\bf s}}_{0,k}^{RRE} }

&fg=000000$;同样分别应用k步Arnoldi和GMRES,从$latex {\left( {{\rm {\bf

I}}-{\rm {\bf A}}} \right){\rm {\bf x}}={\rm {\bf b}}}&fg=000000$中生成

$latex {{\rm {\bf s}}_k^{\mbox{Arnoldi}} }&fg=000000$和$latex {{\rm {\bf

s}}_k^{\mbox{GMRES}} }&fg=000000$。则有如下关系:

$latex {{\rm {\bf s}}_{0,k}^{MPE} ={\rm {\bf s}}_k^{\mbox{Arnoldi}} }

&fg=000000$,$latex {{\rm {\bf s}}_{0,k}^{RRE} ={\rm {\bf s}}_k^{\mbox

{GMRES}} }&fg=000000$。

3.9 循环SMACOF算法

为了加速SMACOF的收敛速度,我们将循环模式应用进来,由于这是个非线性问题,故

外推法的近似极限向量并不一定会有低的stress值。因此我们必须采取某种保护措施



一种方法是检查外推极限的stress值,这个值若较高,则采用最后一次的迭代向量代

替,这个方法在表- 5 循环模式下的向量外推算法中的第五步做一个简单更改即可,

表-6做出了算法描述:
第二种方法是减小步长。

参考文献

[1] Walter .E. Arnoldi. The principle of minimized iterations in the

solution
of the matrix eigenvalue problem. Quart. Appl. Math., 9:17–29, 1951.
[2] Brian T. Bartell, Garrison W. Cottrell, and Richard K. Belew. Latent
semantic indexing is an optimal special case of multidimensional scaling.
In SIGIR ’92: Proceedings of the 15th annual international ACM SIGIR
conference on Research and development in information retrieval, pages
161–167, New York, NY, USA, 1992. ACM Press.
[3] Mikhail Belkin and Partha Niyogi. Laplacian eigenmaps and spectral

tech-
niques for embedding and clustering. In T. G. Dietterich, S. Becker, and
Z. Ghahramani, editors, Advances in Neural Inf. Proc. Sys., volume 14,
pages 585–591, Cambridge, MA, 2002. MIT Press.
[4] Ronald F. Boisvert, Roldan Pozo, Karin Remington, Richard Barrett, and
Jack J. Dongarra. The Matrix Market: A web resource for test matrix
collections. In Ronald F. Boisvert, editor, Quality of Numerical Software,
Assessment and Enhancement, pages 125–137, London, 1997. Chapman
[5] Ingwer Borg and Patrick Groenen. Modern multidimensional scaling: The-
ory and applications. Springer Verlag, New York, 1997.
[6] Ulrik Brandes and Christian Pich. Eigensolver methods for progressive
multidimensional scaling of large data. In Michael Kaufmann and Dorothea
Wagner, editors, Graph Drawing, Karlsruhe, Germany, September 18-20,
2006, pages pp. 42–53. Springer, 2007.
[7] Achi Brandt and Vladimir Mikulinsky. On recombining iterants in multi-
grid algorithms and problems with small islands. SIAM J. Sci. Comput.,
16(1):20–28, 1995.
[8] Alex M. Bronstein, Michael M. Bronstein, and Ron Kimmel. Expression-
invariant face recognition via spherical embedding. In Proc. IEEE Inter-
national Conf. Image Processing (ICIP), 2005.
[9] Alex M. Bronstein, Michael M. Bronstein, and Ron Kimmel. Three-
dimensional face recognition. International Journal of Computer Vision,
64(1):5–30, August 2005.
[10] Michael M. Bronstein, Alex M. Bronstein, R. Kimmel, and I. Yavneh.
Multigrid multidimensional scaling. Numerical Linear Algebra with Appli-
cations, Special issue on multigrid methods, 13(2-3):149–171, March-April
2006.
[11] Stan Cabay and L.W. Jackson. A polynomial extrapolation method for
finding limits and antilimits of vector sequences. SIAM J. Numer. Anal.,
13:734–752, 1976.
[12] Matthew Chalmers. A linear iteration time layout algorithm for

visualising
high-dimensional data. In IEEE Visualization, pages 127–132, 1996.
[13] Ka W. Cheung and Hing C. So. A multidimensional scaling framework for
mobile location using time-of-arrival measurements. IEEE transactions on
signal processing, 53(2):460–470, 2005.
[14] Ronald R. Coifman, Stephane Lafon, Ann B. Lee, Mauro Maggioni, Boaz
Nadler, Frederick Warner, and Steven W. Zucker. Geometric diffusions as
a tool for harmonic analysis and structure definition of data. Proc. Natl.
Acad. Sci. USA, 102(21):7426–7431, May 2005.
[15] Lee G. Cooper. A review of multidimensional scaling in marketing

research.
Applied Psychological Measurement, 7(4):427–450, 1983.
[16] Payel Das, Mark Mol, Hernan Stamati, Lydia E. Kavraki, , and Cecilia
Clementi. Low-dimensional, free-energy landscapes of protein-folding reac-
tions by nonlineardimensionality reduction. Proc. Natl. Acad. Sci. USA,
103(26):9885–9890, June 2006.
[17] Vin de Silva and Joshua B. Tenenbaum. Global versus local methods in
nonlinear dimensionality reduction. In Suzanna Becker, Sebastian Thrun,
and Klaus Obermayer, editors, Advances in Neural Inf. Proc. Sys., pages
705–712. MIT Press, 2002.
[18] Roberto Moreno Diaz and Alexis Quesada Arencibia, editors. Coloring of
DT-MRI Fiber Traces using Laplacian Eigenmaps, Las Palmas de Gran
Canaria, Spain, February 24–28 2003. Springer Verlag.
[19] Robert P. Eddy. Extrapolating to the limit of a vector sequence. In

P.C.C.
Wang, editor, Information Linkage Between Applied Mathematics and In-
dustry, pages 387–396, New York, 1979. Academic Press.
[20] Asi Elad and Ron Kimmel. On bending invariant signatures for surfaces.
IEEE Trans. Pattern Anal. Mach. Intell., 25(10):1285–1295, 2003.
[21] Christos Faloutsos and King-Ip Lin. FastMap: A fast algorithm for

index-
ing, data-mining and visualization of traditional and multimedia datasets.
In Michael J. Carey and Donovan A. Schneider, editors, Proceedings of the
1995 ACM SIGMOD International Conference on Management of Data,
pages 163–174, San Jose, California, 22–25 May 1995.
[22] Ronald A. Fisher. The systematic location of genes by means of

crossover
observations. The American Naturalist, 56:406–411, 1922.
[23] Emden R. Gansner, Yehuda Koren, and Stephen C. North. Graph drawing
by stress majorization. In J´anos Pach, editor, Graph Drawing, volume 3383
of Lecture Notes in Computer Science, pages 239–250. Springer, 2004.
[24] Gene H. Golub and Charles F. Van Loan. Matrix Computations. The Johns
Hopkins University Press, London, third edition, 1996.
[25] Louis Guttman. A general nonmetric technique for finding the smallest
coordinate space for a configuration of points. Psychometrika, 33:469–506,
1968.
[26] Yoshi hiro Taguchi and Yoshitsugu Oono. Relational patterns of gene

ex-
pression via non-metric multidimensional scaling analysis. Bioinformatics,
21(6):730–740(11), march 2005.
[27] Anthony Kearsley, Richard Tapia, and Michael W. Trosset. The solution
of the metric stress and sstress problems in multidimensional scaling using
Newton’s method. Computational Statistics, 13(3):369–396, 1998.
[28] Yosi Keller, Stephane Lafon, and Michael Krauthammer. Protein cluster
analysis via directed diffusion. In The fifth Georgia Tech International
Conference on Bioinformatics, November 2005.
[29] Jospeh B. Kruskal. Multidimensional scaling by optimizing goodness of

fit
to a nonmetric hypothesis. Psychometrika, 29:1–27, 1964.
[30] Cecilio Mar-Molinero and Carlos Serrano-Cinca. Bank failure: a

multidi-
mensional scaling approach. The European Journal of Finance, 7(2):165–
183, June 2001.
[31] M. Me˘sina. Convergence acceleration for the iterative solution of the

equa-
tions X = AX + f. Comput. Methods Appl. Mech. Engrg., 10:165–173,
1977.
[32] Alistair Morrison, Greg Ross, and Matthew Chalmers. Fast multidimen-
sional scaling through sampling, springs and interpolation. Information
Visualization, 2(1):68–77, 2003.
[33] Philip H. Frances Patrick Groenen. Visualizing time-varying

correlations
across stock markets. Journal of Empirical Finance, 7:155–172, 2000.
[34] Robert Pless. Using Isomap to explore video sequences. In Proceedings

of
the 9th International Conference on Computer Vision, pages 1433–1440,
Nice, France, October 2003.
[35] Keith T. Poole. Nonparametric unfolding of binary choice data.

Political
Analysis, 8(3):211–237, March 2000.
[36] Guy Rosman, Alex M. Bronstein, Michael M. Bronstein, and Ron Kimmel.
Topologically constrained isometric embedding. In Proc. Conf. on Machine
Learning and Pattern Recognition (MLPR), 2006.
[37] Sam T. Roweis and Lawrence K. Saul. Nonlinear dimensionality reduction
by locally linear embedding. Science, 290:2323–2326, 2000.
[38] Yousef Saad and Martin H. Schultz. GMRES: A generalized minimal resid-
ual method for solving nonsymmetric linear systems. SIAM J. Sci. Statist.
Comput., 7:856–869, 1986.
[39] J.R. Schmidt. On the numerical solution of linear simultaneous

equations
by an iterative method. Phil. Mag., 7:369–383, 1941.
[40] Eric. L. Schwartz, Alan Shaw, and Estarose Wolfson. A numerical

solution
to the generalized mapmaker’s problem: Flattening nonconvex polyhedral
surfaces. IEEE Trans. Pattern Anal. Mach. Intell., 11:1005–1008, Novem-
ber 1989.
[41] Daniel Shanks. Nonlinear transformations of divergent and slowly

conver-
gent sequences. J. Math. and Phys., 34:1–42, 1955.
[42] Avram Sidi. Convergence and stability properties of minimal polyno-
mial and reduced rank extrapolation algorithms. SIAM J. Numer. Anal.,
23:197–209, 1986. Originally appeared as NASA TM-83443 (1983).
[43] Avram Sidi. Extrapolation vs. projection methods for linear systems of
equations. J. Comp. Appl. Math., 22:71–88, 1988.
[44] Avram Sidi. Efficient implementation of minimal polynomial and reduced
rank extrapolation methods. J. Comp. Appl. Math., 36:305–337, 1991.
Originally appeared as NASA TM-103240 ICOMP-90-20.
[45] Avram Sidi. Convergence of intermediate rows of minimal polynomial and
reduced rank extrapolation tables. Numer. Algorithms, 6:229–244, 1994.
[46] Avram Sidi and J. Bridger. Convergence and stability analyses for some
vector extrapolation methods in the presence of defective iteration

matrices.
J. Comp. Appl. Math., 22:35–61, 1988.
[47] Avram Sidi, William F. Ford, and David A. Smith. Acceleration of con-
vergence of vector sequences. SIAM J. Numer. Anal., 23:178–196, 1986.
Originally appeared as NASA TP-2193, (1983).
[48] Avram Sidi and Yair Shapira. Upper bounds for convergence rates of

vector
extrapolation methods on linear systems with initial iterations. Technical
Report 701, Computer Science Department, Technion–Israel Institute of
Technology, 1991. Appeared also as NASA Technical memorandum 105608,
ICOMP-92-09, (1992).
[49] Avram Sidi and Yair Shapira. Upper bounds for convergence rates of ac-
celeration methods with initial iterations. Numer. Algorithms, 18:113–132,
1998.
[50] Stig Skelboe. Computation of the periodic steady-state response of

nonlin-
ear networks by extrapolation methods. IEEE Trans. Circuits and Systems,
27:161–175, 1980.
[51] David A. Smith, William F. Ford, and Avram Sidi. Extrapolation methods
for vector sequences. SIAM Rev., 29:199–233, 1987.
[52] Joshua B. Tenenbaum, Vin de Silva, and John C. Langford. A global
geometric framework for nonlinear dimensionality reduction. Science,
290(5500):2319–2323, December 2000.
[53] Warren S. Torgerson. Multidimensional scaling I. theory and method.
PSym, 17:401–419, 1952.
[54] Shusaku Tsumoto and Shoji Hirano. Visualization of rule’s similarity

using
multidimensional scaling. In ICDM, pages 339–346, 2003.
[55] Jason Tsong-Li Wang, Xiong Wang, King-Ip Lin, Dennis Shasha, Bruce A.
Shapiro, and Kaizhong Zhang. Evaluating a class of distance-mapping
algorithms for data mining and clustering. In KDD ’99: Proceedings of the
fifth ACM SIGKDD international conference on Knowledge discovery and
data mining, pages 307–311, New York, NY, USA, 1999. ACM Press.
[56] Kilian Q. Weinberger and Laurence K. Saul. Unsupervised learning of
image manifolds by semidefinite programming. In Proceedings of the IEEE
Conference on Computer Vision and Pattern Recognition, volume 2, pages
988–995, Washington D.C., 2004. IEEE Computer Society.
[57] Tynia Yang, Jinze Liu, Leonard McMillan, and Wei Wang. A fast approx-
imation to multidimensionalscaling. In IEEE workshop on Computation
Intensive Methods for Computer Vision, 2006.
[58] Gil Zigelman, Ron Kimmel, and Nahum Kiryati. Texture mapping us-
ing surface flattening via multidimensional scaling. IEEE Transactions on
Visualization and Computer Graphics, 8(2):198–207, 2002.