一、基本技術

1)MATLAB索引或引用(MATLAB Indexing or Referencing)

在MATLAB中有三種基本方法可以選取一個矩陣的子陣。它們分別是下標法，線性法和邏輯法(subscripted, linear,

and

logical)。

1.1)下標法

非常簡單，看幾個例子就好。

A = 6:12;

A([3,5])

ans =

8 10

A([3:2:end])

ans =

8 10 12

A = [11 14 17;12 15 18;13 16

19];

A(2:3,2)

ans =

15

16

1.2)線性法

二維矩陣以列優(yōu)先順序可以線性展開，可以通過現(xiàn)行展開后的元素序號來訪問元素。

A = [11 14 17;12 15 18;13 16

19];

A(6)

ans =

16

A([3,1,8])

ans =

13 11 18

A([3;1;8])

ans =

13

11

18

1.3)邏輯法

用一個和原矩陣具有相同尺寸的0-1矩陣，可以索引元素。在某個位置上為1表示選取元素，否則不選。得到的結果是一個向量。

A = 6:10;

A(logical([0 0 1 0 1]))

ans =

8 10

A =[1 2;3 4];

B = [1 0 0 1];

A(logical(B))

ans =

1?4

2)數(shù)組操作和矩陣操作(Array Operations vs. Matrix

Operations)

對矩陣的元素一個一個孤立進行的操作稱作數(shù)組操作；而把矩陣視為一個整體進行的運算則成為矩陣操作。MATLAB運算符*,/,\,^都是矩陣運算，而相應的數(shù)組操作則是.*,

./, .\, .^

A=[1 0 ;0 1];

B=[0 1 ;1 0];

A*B % 矩陣乘法

ans =

0?1

1?0

A.*B?% A和B對應項相乘

ans =

0?0

0?0

3)布朗數(shù)組操作(Boolean Array

Operations)

對矩陣的比較運算是數(shù)組操作，也就是說，是對每個元素孤立進行的因此其結果就不是一個“真”或者“假”，而是一堆“真假”。這個結果就是布朗數(shù)組。

D = [-0.2 1.0 1.5 3.0 -1.0 4.2

3.14];

D >= 0

ans =

0 1 1 1 0 1

1

如果想選出D中的正元素：

D = D(D>0)

D =

1.0000?1.5000?3.0000?4.2000?3.1400

除此之外，MATLAB運算中會出現(xiàn)NaN,Inf,-Inf。對它們的比較參見下例

Inf==Inf返回真

Inf<1返回假

NaN==NaN返回假

同時，可以用isinf，isnan判斷，用法可以顧名思義。

在比較兩個矩陣大小時，矩陣必須具有相同的尺寸，否則會報錯。這是你用的上size和isequal,isequalwithequalnans(R13及以后)。

4)從向量構建矩陣(Constructing Matrices from

Vectors)

在MATLAB中創(chuàng)建常數(shù)矩陣非常簡單，大家經(jīng)常使用的是：

A =

ones(5,5)*10

但你是否知道，這個乘法是不必要的？

A = 10;

A = A(ones(5,5))

A =

10 10 10 10 10

10 10 10 10 10

10 10 10 10 10

10 10 10 10 10

10 10 10 10

10

類似的例子還有：

v = (1:5)';

n = 3;

M = v(:,ones(n,1))

M =

1 1 1

2 2 2

3 3 3

4 4 4

5 5 5

事實上，上述過程還有一種更加容易理解的實現(xiàn)方法：

A = repmat(10,[5 5]);

M = repmat([1:5]', [1,3]);

其中repmat的含義是把一個矩陣重復平鋪，生成較大矩陣。

更多詳細情況，參見函數(shù)repmat和meshgrid。

5)相關函數(shù)列表(Utility Functions)

ones 全1矩陣

zeros?全0矩陣

reshape?修改矩陣形狀

repmat?矩陣平鋪

meshgrid?3維plot需要用到的X－Y網(wǎng)格矩陣

ndgrid?n維plot需要用到的X－Y－Z...網(wǎng)格矩陣

filter?一維數(shù)字濾波器，當數(shù)組元素前后相關時特別有用。

cumsum?數(shù)組元素的逐步累計

cumprod?數(shù)組元素的逐步累計

eye?單位矩陣

diag?生成對角矩陣或者求矩陣對角線

spdiags?稀疏對角矩陣

gallery?不同類型矩陣庫

pascal?Pascal

矩陣

hankel?Hankel

矩陣

toeplitz?Toeplitz

矩陣

二、擴充的例子

6)作用于兩個向量的矩陣函數(shù)

假設我們要計算兩個變量的函數(shù)F

F(x,y) = x*exp(-x^2 -

y^2)

我們有一系列x值，保存在x向量中，同時我們還有一系列y值。

我們要對向量x上的每個點和向量y上的每個點計算F值。換句話說，我們要計算對于給定向量x和y的所確定的網(wǎng)格上的F值。

使用meshgrid，我們可以復制x和y來建立合適的輸入向量。然后可以使用第2節(jié)中的方法來計算這個函數(shù)。

x = (-2:.2:2);

y = (-1.5:.2:1.5)';

[X,Y] = meshgrid(x, y);

F = X .* exp(-X.^2 -

Y.^2);

如果函數(shù)F具有某些性質，你甚至可以不用meshgrid，比如

F(x,y) = x*y ，則可以直接用向量外積

x = (-2:2);

y = (-1.5:.5:1.5);

x'*y

在用兩個向量建立矩陣時，在有些情況下，稀疏矩陣可以更加有效地利用存儲空間，并實現(xiàn)有效的算法。我們將在第8節(jié)中以一個實例來進行更詳細地討論.

7)排序、設置和計數(shù)(Ordering, Setting, and Counting

Operations)?在迄今為止討論過的例子中，對向量中一個元素的計算都是獨立?于同一向量的其他元素的。但是，在許多應用中，你要做的計算?則可能與其它元素密切相關。例如，假設你用一個向量x來表示一?個集合。不觀察向量的其他元素，你并不知道某個元素是不是一?個冗余元素，并應該被去掉。如何在不使用循環(huán)語句的情況下刪除?冗余元素，至少在現(xiàn)在，并不是一個明顯可以解決的問題。

解決這類問題需要相當?shù)闹乔伞Ｒ韵陆榻B一些可用的基本工具

max 最大元素

min 最小元素

sort?遞增排序

unique?尋找集合中互異元素(去掉相同元素)

diff?差分運算符[X(2) - X(1), X(3) - X(2), ... X(n) -

X(n-1)]

find?查找非零、非NaN元素的索引值

union?集合并

intersect?集合交

setdiff?集合差

setxor?集合異或

繼續(xù)我們的實例，消除向量中的多余元素。注意：一旦向量排序后，?任何多余的元素就是相鄰的了。同時，在任何相等的相鄰元素在向量?diff運算時變?yōu)榱恪＿@是我們能夠應用以下策略達到目的。我們現(xiàn)在?在已排序向量中，選取那些差分非零的元素。

% 初次嘗試。不太正確!

x = sort(x(:));

difference = diff(x);

y = x(difference~=0);

這離正確結果很近了，但是我們忘了diff函數(shù)返回向量的元素個數(shù)比?輸入向量少1。在我們的初次嘗試中，沒有考慮到最后一個元素也可能?是相異的。為了解決這個問題，我們可以在進行差分之前給向量x加入?一個元素，并且使得它與以前的元素一定不同。一種實現(xiàn)的方法是增?加一個NaN。

% 最終的版本。

x = sort(x(:));

difference = diff([x;NaN]);

y =

x(difference~=0);

我們使用(:)運算來保證x是一個向量。我們使用~=0運算，而不用find?函數(shù)，因為find函數(shù)不返回NaN元素的索引值，而我們操作中差分的最?后元素一定是NaN。這一實例還有另一種實現(xiàn)方式：

y=unique(x);

后者當然很簡單，但是前者作為一個練習并非無用，它是為了練習使用?矢量化技術，并示范如何編寫你自己的高效代碼。此外，前者還有一個?作用：Unique函數(shù)提供了一些超出我們要求的額外功能，這可能降低代?碼的執(zhí)行速度。

假設我們不只是要返回集合x，而且要知道在原始的矩陣里每個相異元素?出現(xiàn)了多少個“復本”。一旦我們對x排序并進行了差分，我們可以用?find來確定差分變化的位置。再將這個變化位置進行差分，就可以得到?復本的數(shù)目。這就是"diff

of find of

diff"的技巧。基于以上的討論，?我們有：

% Find the redundancy in a vector

x

x = sort(x(:));

difference =

diff([x;max(x)+1]);

count =

diff(find([1;difference]));

y = x(find(difference));

plot(y,count)

這個圖畫出了x中每個相異元素出現(xiàn)的復本數(shù)。注意，在這里我們避開了?NaN，因為find不返回NaN元素的索引值。但是，作為特例，NaN和Inf?的復本數(shù)可以容易地計算出來：

count_nans =

sum(isnan(x(:)));

count_infs =

sum(isinf(x(:)));

另一個用于求和或者計數(shù)運算的矢量化技巧是用類似建立稀疏矩陣的方?法實現(xiàn)的。這還將在第9節(jié)中作更加詳細的討論.

8)稀疏矩陣結構(Sparse Matrix

Structures)

在某些情況下，你可以使用稀疏矩陣來增加計算的效率。如果你構造一?個大的中間矩陣，通常矢量化更加容易。在某些情況下，你可以充分利?用稀疏矩陣結構來矢量化代碼，而對于這個中間矩陣不需要大的存儲空?間。

假設在上一個例子中，你事先知道集合y的域是整數(shù)的子集，

{k+1,k+2,...k+n}；即，

y = (1:n) +

k

例如，這樣的數(shù)據(jù)可能代表一個調色板的索引值。然后，你就可以對集?合中每個元素的出現(xiàn)進行計數(shù)(構建色彩直方圖？譯者)。這是對上一?節(jié)中"diff

of find of diff"技巧的一種變形。

現(xiàn)在讓我們來構造一個大的m x

n矩陣A，這里m是原始x向量中的元素數(shù)，?n是集合y中的元素數(shù)。

A(i,j) = 1 if x(i) =

y(j)

0 otherwise

回想一下第3節(jié)和第4節(jié)，你可能認為我們需要從x和y來構造矩陣A。如果?當然可以，但要消耗許多存儲空間。我們可以做得更好，因為我們知道，?矩陣A中的多數(shù)元素為0，x中的每個元素對應的行上只有一個值為1。

以下就是構造矩陣的方法(注意到y(tǒng)(j) =

k+j，根據(jù)以上的公式)：

x = sort(x(:));

A = sparse(1:length(x), x+k, 1, length(x),

n);

現(xiàn)在我們對A的列進行求和，得到出現(xiàn)次數(shù)。

count =

sum(A);

在這種情況下，我們不必明確地形成排序向量y，因為我們事先知道

y = 1:n + k.

這里的關鍵是使用數(shù)據(jù)，(也就是說，用x控制矩陣A的結構)。由于x在?一個已知范圍內取整數(shù)值，我們可以更加有效地構造矩陣。

假設你要給一個很大矩陣的每一列乘以相同的向量。使用稀疏矩陣，不僅?可以節(jié)省空間，并且要比在第5節(jié)介紹的方法更加快速.

下面是它的工作?方式：

F = rand(1024,1024);

x = rand(1024,1);

% 對F的所有行進行點型乘法.

Y = F * diag(sparse(x));

% 對F的所有列進行點型乘法.

Y = diag(sparse(x)) *

F;

我們充分利用矩陣乘法算符來執(zhí)行大規(guī)模運算，并使用稀疏矩陣以防止臨?時變量變得太大。

9)附加的例子(Additional

Examples)

下面的例子使用一些在本技術手冊中討論過的方法，以及其它一些相關方?法。請嘗試使用tic

和toc(或t=cputime和cputime-t)，看一下速度加快?的效果。

用于計算數(shù)組的雙重for循環(huán)。

使用的工具：數(shù)組乘法

優(yōu)化前：

A = magic(100);

B = pascal(100);

for j =

1:100

for k =

1:100;

X(j,k) = sqrt(A(j,k)) * (B(j,k) -

1);

end

end

優(yōu)化后：

A = magic(100);

B = pascal(100);

X = sqrt(A).*(B-1);

用一個循環(huán)建立一個向量，其元素依賴于前一個元素

使用的工具：FILTER, CUMSUM,

CUMPROD

優(yōu)化前：

A = 1;

L = 1000;

for i = 1:L

A(i+1) = 2*A(i)+1;

end

優(yōu)化后：

L = 1000;

A = filter([1],[1

-2],ones(1,L+1));

如果你的向量構造只使用加法或乘法，你可使用cumsum或cumprod函數(shù)。

優(yōu)化前：

n=10000;

V_B=100*ones(1,n);

V_B2=100*ones(1,n);

ScaleFactor=rand(1,n-1);

for i =

2:n

V_B(i) =

V_B(i-1)*(1+ScaleFactor(i-1));

end

for

i=2:n

V_B2(i) = V_B2(i-1)+3;

end

優(yōu)化后：

n=10000;

V_A=100*ones(1,n);

V_A2 = 100*ones(1,n);

ScaleFactor=rand(1,n-1);

V_A=cumprod([100

1+ScaleFactor]);

V_A2=cumsum([100

3*ones(1,n-1)]);

向量累加，每5個元素進行一次：

工具：CUMSUM , 向量索引

優(yōu)化前：

% Use an arbitrary vector,

x

x = 1:10000;

y = [];

for n =

5:5:length(x)

y = [y sum(x(1:n))];

end

優(yōu)化后(使用預分配)：

x = 1:10000;

ylength = (length(x) -

mod(length(x),5))/5;

% Avoid using ZEROS command during

preallocation

y(1:ylength) = 0;

for n =

5:5:length(x)

y(n/5) = sum(x(1:n));

end

優(yōu)化后(使用矢量化,不再需要預分配)：

x = 1:10000;

cums = cumsum(x);

y = cums(5:5:length(x));

操作一個向量，當某個元素的后繼元素為0時，重復這個元素：

工具：FIND, CUMSUM,

DIFF

任務：我們要操作一個由非零數(shù)值和零組成的向量，要求把零替換成為?它前面的非零數(shù)值。例如，我們要轉換下面的向量：

a=2; b=3; c=5; d=15; e=11;

x = [a 0 0 0 b 0 0 c 0 0 0 0 d 0 e 0 0 0 0

0];

為：

x = [a a a a b b b c c c c c d d e e e e e

e];

解(diff和cumsum是反函數(shù))：

valind = find(x);

x(valind(2:end)) =

diff(x(valind));

x =

cumsum(x);

將向量的元素累加到特定位置上

工具：SPARSE

優(yōu)化前：

% The values we are summing at designated

indices

values = [20 15 45 50 75 10 15 15 35 40

10];

% The indices associated with the values are summed

cumulatively

indices = [2 4 4 1 3 4 2 1 3 3

1];

totals =

zeros(max(indices),1);

for n =

1:length(indices)

totals(indices(n)) = totals(indices(n)) +

values(n);

end

優(yōu)化后：

indices = [2 4 4 1 3 4 2 1 3 3

1];

totals =

full(sparse(indices,1,values));

注意：這一方法開辟了稀疏矩陣的新用途。在使用sparse命令創(chuàng)建稀疏矩陣?時，它是對分配到同一個索引的所有值求和，而不是替代已有的數(shù)值。這稱?為"向量累加"，是MATLAB處理稀疏矩陣的方式。

總結

以上是生活随笔為你收集整理的matlab 矢量化,matlab矢量化编程简要的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內容還不錯，歡迎將生活随笔推薦給好友。