CodyNote016：「拉鏈矩陣」的構造

來自專欄 Cody習題6 人贊了文章

馬良，祁彬彬

序言

儘管有「矩陣實驗室」的美名，但MATLAB操作矩陣對角線元素的能力還是個軟肋和痛點。當然，這是相對於MATLAB本身，其他矩陣元素特徵的操縱能力而言。本篇打算探討矩陣元素的「拉鏈」式提取，這仍然是綜合基本矩陣元素操作、索引控制以及特殊矩陣形式構造的綜合問題。

可能有人不以為然，認為既然MATLAB內部已經有diag、spdiags、blkdiag、eye等一眾犀利武器，對角線矩陣元素的操作那還不是360度無死角，上下左右隨意玩弄？事實上考慮到一些攪屎棍級的障礙，例如元素提取順序，這類問題也並非想像般輕鬆地手到擒來。

「拉鏈矩陣」系列問題

Pro1167. matrix zigzag

【鏈接】：https://ww2.mathworks.cn/matlabcentral/cody/problems/1167

【原題】：Unfold a 2-D matrix to a 1-D array in zig-zag order, e.g., for matrix

[ 1 2 3 ; 4 5 6 ; 7 8 9 ]

the resulting 1-D array should be

[ 1 2 4 7 5 3 6 8 9 ]

【解釋】：前面提到對角線元素，加上所謂「拉鏈」矩陣，意指沿對角線或者反對角線方向，左一下、右一下的走「之」字提取元素，本例是從輸入矩陣第1個元素開始，沿反對角線一正一反，「之」字形提取原矩陣元素，並把提取元素轉換為元素總數相同的1-D數組。

Test

Code Input

1

%%x = [1 2; 3 4];y_correct = [1 2 3 4];assert(isequal(zigzag(x),y_correct))

2

%%x = [ 1 2 3; 4 5 6; 7 8 9];y_correct = [ 1 2 4 7 5 3 6 8 9];assert(isequal(zigzag(x),y_correct))

3

%%x = magic(4);y_correct = [16 2 5 9 11 3 13 10 7 4 14 6 8 12 15 1];assert(isequal(zigzag(x),y_correct))

【分析】：照原題例子，按反對角線，一正序一倒序提取元素，結果依序為：

[[1],[2 4],[7 5 3],[6 8],[9]]

Test Suite也不複雜，不需要搞「空矩陣」之類的設定，難度僅體現在「對角線」上。

逐個元素提取的for循環+if判定

如果逐個元素提取，那麼大量判斷、行列索引號變換+1或-1是很容易預料到的，例如：

function y = zigzag(x) n = size(x,1); i = 1; j = 1; k = 1; y = reshape(x,1,[]); d = true; % true for up-right, false for down-left while i<n | j<n k = k+1; if d % up right if i==1 | j==n % hit top row or right column, now move down-left d = false; if j==n i = i+1; else j = j+1; end else % continue up-right i = i - 1; j = j + 1; end elseif i==n | j==1 % hit bottom row or left column, now move up-right d = true; if i==n j = j+1; else i = i+1; end else % continue down-left i = i+1; j = j-1; end y(k) = x(i,j); endend

用到MATLAB函數很少，除控制流程語句for、while、if等，只有size和reshape，多數情況下這意味著沒有全部發揮出MATLAB自身的潛力。

逐條對角線元素提取

比逐個元素提取更進一步的是按對角線提取，因為對角線類型問題，容易想到diag命令，diag不僅能提取主對角線元素，也可由改變第2參數大小提取其他對角線；同時，題意要求反對角線，rot90、flip等命令即成為輔助函數備選。比較典型的是如下這種：

function y = zigzag(x) x = rot90(x); r = size(x,1); y = []; flag = true; for i = -r + 1 : r - 1 if flag y = [y,fliplr(diag(x,i))]; else y = [y,diag(x,i)]; end flag = not(flag); endend

利用循環來動態改變flag數值，以操縱輸入矩陣中提取的元素順序。多數初學者會選擇接受使用此演算法寫出求解代碼。上述演算法還有一定優化空間，比如我略修改Tim另一問題中相同思路，提交了一個方案：

function ans = zigzag(a) a = rot90(a); []; for k = 1:2*size(a,2)-1 [ans flip(diag(a,k-size(a,1)), 1+mod(k,2))]; end end

以上2個方案基本屬於同樣思路，但Tim的構造方式更緊湊，用mod對2求余動態獲得對角線元素翻轉控制的flag數值，替換大段if語句，提取對角線第2參數則用「k-size(a,1)」一步到位，這是個堪稱精巧的寫法，用簡單四則運算構造替換if語句，這是高手常用的縮減size手法，值得借鑒。此外還有劉鵬把rot90寫進循環的辦法：

function ans = zigzag(x) n = length(x); ; for m = 1:2*n-1 [ans rot90(diag(rot90(x),m-n)., 2*mod(m,2))]; endend

循環節「1:2*size(a,2)-1」代表矩陣對角線條數，但和diag第2參數有差異，以x=magic(4)為例，m或k=1時，diag第2參數數值應該用m-size(x,1)=1-4=-3，這有些麻煩，因此還有一種方式是在循環節上直接改成第2參數：

function ans=zigzag(x)[]; n=length(x)-1;for i=n:-1:-n x=x; [ans; diag(fliplr(x),i)];end

Khaled Hamed用「i=n:-1:-n」避免diag第2參數在循環內的運算，在循環中轉置和左右翻轉整個矩陣獲得順序和逆序對角線元素，即：變化的不是對角線，是整個矩陣。這和劉鵬以及之前Tim的思路均有部分重合之處。而J.S.Kowontan提交代碼則類似劉鵬和Khaled Hamed兩種計算方法的綜合：

function ans = zigzag(x) length(x)-1; cell2mat(arrayfun(@(z)rot90(diag(rot90(x),z),2*(z+ans-1)),-ans:ans,uni,0));end

當然，也可以搞成「一行流」，這是無所謂，看個人喜好了：

function ans = zigzag(x)cell2mat(arrayfun(@(ii)rot90(diag(fliplr(x),-ii),2*(ii+length(x))),-length(x):length(x),un,0));end

最後，Khaled Hamed的演算法需要輸入x及動態改變長度的輸出數組不停轉置，這是因為diag輸出為列向量的原因，無論如何，還是有點逼死強迫症，所以Bryant Tran用rot90替換轉置：

function ans = zigzag(x) []; for ii = -length(x):length(x) [rot90(diag(fliplr(x),ii),2*(ii+length(x))+1) ans]; endend

以上代碼方案思路基本相同，都靠diag提取對角線元素，再翻轉矩陣或對角線獲得元素的正確順序，翻轉對角線元素時，依據mod對2求余是個非常棒的主意，一個命令就去掉了大段的if判斷。

稀疏矩陣提取對角線元素的spdiags命令

本節的思路和CodyNote015提到的spdiags有關，這是用於稀疏矩陣中提取存在非零元素對角線元素的命令。實際上這也是本問題中縮小size的最有效手段，因為spdiags能把不好進行索引操控的對角線變成列元素，關於spdiags也可以參看我和彬彬書里關於連續質數對角線的例子。

function ans = zigzag(x) spdiags(flipud(x)); ans(:,2:2:end)=flipud(ans(:,2:2:end)); nonzeros(ans);end

以x=magic(4)為例，代碼第1句把所有對角線元素變成列向量，問題轉化為列元素的提取：

>> spdiags(flipud(x))ans = 16 5 9 4 0 0 0 0 2 11 7 14 0 0 0 0 3 10 6 15 0 0 0 0 13 8 12 1

下一句把需要逆序提取的元素用flipud翻轉，最後nonzeros一次提取非零元素，此時元素順序就是正確的。

這個演算法也不算完美，因為x本身有0元素時，這段代碼執行時會出現問題。只是Test Suite沒出現存在0元素的例子而逃過一劫。

最後從本例代碼出發深入一步，談談另一個有趣的問題：R2016b版本後，隱式擴展操作的問題。

由前述，可用mod命令替代if語句，保證對角線元素提取時，「正、逆、正...」的順序，剛才Alfonso代碼中用到翻轉2:2:end的索引操控方法保證元素順序，那麼，能不能結合mod和spdiags寫一個代碼呢？不妨試一試：

>> a=magic(4)a = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1>> spdiags(rot90(a))ans = 16 2 3 13 0 0 0 0 5 11 10 8 0 0 0 0 9 7 6 12 0 0 0 0 4 14 15 1>> rot90(ans,2*mod(1:2*length(a)-1,2))Error using rot90 (line 23)k must be a scalar.

我們發現似乎無法實現原來的意圖，因為rot90的第2參數k是90°的整數倍，且只接受單一數值(scalar)，而「2*mod(1:2*length(a)-1,2)」的輸出為數組：

>> 2*mod(1:2*length(a)-1,2)ans = 2 0 2 0 2 0 2

因此運行提示錯誤。所以...真的不能用spdiags+mod實現問題的代碼嗎？

回答問題前，再秀一把bsxfun，這是R2016b版本之前完成隱式擴展的函數，用法靈活且執行效率極高，一眾高手對其相當推崇，但R2016b推出後，常規矩陣操作中合併了隱式擴展功能，似乎二者間就不再作區分，因此這讓很多人...好吧，主要是我自己，感覺原版本bsxfun函數似乎已可以徹底退出歷史舞台了，但這個問題的解決代碼證明這是毫無根據的錯覺：

function ans = zigzag(a)nonzeros(bsxfun(@rot90,spdiags(rot90(a)),2*mod(1:2*length(a)-1,2)));end

由於句柄「@rot90」調用第2參數為1×7數組，但注意第1參數spdiags(rot90(a))運行結果為4×7矩陣，二者列數相同，因此自動執行隱式擴展，每列完成rot90操作，k值在奇數列為2，對列向量而言，這就是上下翻轉，偶數列為0則保持元素位置不變，以調用句柄rot90的方式，乾淨利索完成對奇數列的翻轉操作。

這是一段極度具有啟發性的代碼，提示了bsxfun函數仍具有普通矩陣操作命令所不具備的功能，也就是深入到函數內部對參數結構做隱式擴展。

利用特殊矩陣hilb的擴展構造

LY Cao給出了一個不落俗套，甚至相當妖艷的方案：通過特殊矩陣hilb做底部構造，再由試湊獲得結果。

function ans = zigzag(x) 1./hilb(length(x)); x(sortrowsc(nonzeros(ans+(-1).^ans./sum(ans)),1));end

有2點值得注意：首先使用了特殊矩陣構造輸出變數的大體輪廓，以t=magic(6)為例：

>> t=magic(6)t = 35 1 6 26 19 24 3 32 7 21 23 25 31 9 2 22 27 20 8 28 33 17 10 15 30 5 34 12 14 16 4 36 29 13 18 11>> 1./hilb(length(t))ans = 1 2 3 4 5 6 2 3 4 5 6 7 3 4 5 6 7 8 4 5 6 7 8 9 5 6 7 8 9 10 6 7 8 9 10 11

這實際上為形成索引矩陣提供了備選，如遇到類似矩陣元素特徵，可參考類似的構造思路。

其次，代碼使用help中無法查到的sortrowsc命令，。通俗地講，這是有著M-Function外包裝的C-MEX函數，它是sortrows的輔助函數，其結果返回sortrows的索引，但輸入必須是列向量：

>> a=[16 5 9 4]a = 16 5 9 4>> sortrows(a,1)ans = 4 5 9 16>> sortrowsc(a,1)ans = 4 2 3 1

在這段代碼中，顯然並不出於效率目的，而是為直接輸出sortrows索引，如採用sortrows或者sort，則需要用如下代碼增加一行：

>> [~,idx]=sortrows(magic(4),1)idx = 4 2 3 1

由於代碼中sortrowsc返回的是索引，因此對於矩陣存在非零元素的情況，也同樣適用，這點不同於spdiags+nonzeros+rot90的方案。

Pro44102.Stop ZigZag scanning N*N Matrix at any diag you want

【鏈接】：https://ww2.mathworks.cn/matlabcentral/cody/problems/44102

【題目】：Suppose that we have a 2-D matrix and we try to obtain a 1-D array in zig-zag order, but not all values of our 2-D matrix e.g., for matrix

x= [1 2 3; 4 5 6; 7 8 9]

the resulting 1-D array should be:

stop_zig(x,3)=[ 1 2 4 7 5 3]stop_zig(x,1)=1stop_zig(x,-1)=9stop_zig(x,-2)=[9 8 6]

【解釋】：仍然是拉鏈矩陣元素索引，不過改成提取部分元素，且提取時的方向根據第2參數n的取值是否大於0，存在正向和逆向倒序兩種。例如n=3，代表從左上角開始，提取前3條對角線；如果n=-3，代表從右下角起，自右下至左上，提取3條對角線。

Test

Code Input

1

%%x = 1; n=1;y_correct = 1;assert(isequal(stop_zig(x,n),y_correct))

2

%%x =[1 2 3; 4 5 6; 7 8 9]; n=3;y_correct = [ 1 2 4 7 5 3];assert(isequal(stop_zig(x,n),y_correct))

3

%%x =[1 2 3; 4 5 6; 7 8 9]; n=-3;y_correct = [ 9 8 6 3 5 7];assert(isequal(stop_zig(x,n),y_correct))

4

%%x =[1 2 3; 4 5 6; 7 8 9]; n=2;y_correct = [ 1 2 4];assert(isequal(stop_zig(x,n),y_correct))

5

%%x =[1 2 3 4; 4 5 6 7; 7 8 9 10; 10 11 12 13]; n=-3;y_correct = [ 13 12 10 7 9 11];assert(isequal(stop_zig(x,n),y_correct))

【分析】：兩個題目思路類似，spdiags或者mod等等都可以用，只是要知道兩件事：沿什麼方向從哪兒開始，以及何時停止？

旋轉、對角線提取和省掉的if判斷

Test Suite全為非0元素，因此採用spdiags+rot90組合的方案相對方便：

function ans = stop_zig(x,n) spdiags(rot90(x,1)); if n<0 fliplr(ans); ans(:,2:2:end)=flipud(ans(:,2:2:end)); else ans(:,1:2:end)=flipud(ans(:,1:2:end)); end nonzeros(ans(:,1:abs(n))); end

先用spdiags+rot90提取對角線元素，構造主體框架，由於翻轉列的位置可能從x(1)或x(end)開始，所以對矩陣整體依據第2輸入n>0或n<0做翻轉與否的判斷，因此引起反轉列起始位置也有所不同，判斷執行語句中相應有奇偶列之分。總體上講所需函數都提到了，但寫得還是略啰嗦，體現在2個方面:if語句和兩個判定分支中大多數重複的翻轉語句。怎麼改呢？按前述，if語句可通過邏輯數組替換，這裡顯然可由「n<0」做點兒文章：

function ans = stop_zig(x,n)flip(spdiags(rot90(x,1)),3-(n<0));t=(n<0)+1:2:size(ans,2);ans(:,t)=flipud(ans(:,t));nonzeros(ans(:,1:abs(n))); end

依輸入n的不同，「(n<0)」可能是0或1，變數t中的(n<0)+1相應為1或2，原來的翻轉語句合併為1個，在n<0時左右翻轉矩陣「spdiags(rot90(x,1))」，flip中通過「3-(n<0)」獲得不同情況下的翻轉。

彬彬採用整體處理的思路：

function ans = stop_zig(x,n)spdiags(flip(x))ans(:,2:2:end)=flipud(ans(:,2:2:end));nonzeros(ans);if n < 0 flip(ans);endans(1:nchoosek(abs(n)+1,2));end

代碼中翻轉的都是偶數列，這部分未做判斷，但得到的非零元素行，如果第2輸入參數n<0，此時再左右翻轉，最後提取的總是前n條對角線元素。

順便說一句，因為對角線元素是1-2-3-...的規律，所以提取元素的數量就是1-3-6-...，這正好是pascal三角形的規律，nchoosek正好派上用場。

多次翻轉與索引操控

Alfonso通過常用函數，規規矩矩的索引操控和構造解決了這個問題，雖然命令常見，還是一如既往燒腦。

function y = stop_zig(x,n) if n<0, x=flipud(fliplr(x)); end [i,j]=find(1|x); k=i+j; [~,idx]=sort(k+eps*j.*(-1).^k); y=x(idx(k(idx)<=abs(n)+1));end

本質上這段代碼和前面LY Cao的計算方法相同，分析如下：

(1). 第1行採用改變整個矩陣配合自矩陣首位搜索的順序，不過if判定可通過構造rot90第2參數的邏輯運算省略，也就是把

if n<0 x=flipud(fliplr(x)); end

換成：

x=rot90(x,2*(n<0))

(2). 後面兩句所得結果等效於LY Cao代碼中用hilb構造的矩陣變成a(:)的形式，即：

>> xx = 1 2 3 4 5 6 7 8 9>> [i,j]=find(1|x)i = 1 2 3 1 2 3 1 2 3j = 1 1 1 2 2 2 3 3 3>> k=i+jk = 2 3 4 3 4 5 4 5 6

結果等效於：

>> k1=1./hilb(3)+1k1 = 2 3 4 3 4 5 4 5 6>> k1(:)ans = 2 3 4 3 4 5 4 5 6

題外話：如果單純只需要獲得類似上述k1變數的形式，更簡單辦法是隱式擴展：

>> [1 2 3]+[1 2 3]ans = 2 3 4 3 4 5 4 5 6

或者bsxfun的版本：

>> bsxfun(@plus,1:3,(1:3))ans = 2 3 4 3 4 5 4 5 6

但find輸出參數i和j(行列索引號)在後續試湊代碼中還要用，硬套隱式擴展就繞圈兒了。

(3). 索引排序部分，即：「[~,idx]=sort(k+eps*j.*(-1).^k)」通過試湊構造矩陣，其升序排列的元素在原矩陣的索引idx等效於拉鏈矩陣的檢索順序，eps顯然是為決定相等的行列索引的大小，從拉鏈矩陣對角線元素提取到這樣的試湊矩陣索引，真的是很難想到的內在數據關係。

逐條對角線提取翻轉的循環方式

yurenchu的代碼方案如下：

function ans = stop_zig(x,n) x = rot90(x,sign(n)); []; for k = 1:abs(n) [ans flip(diag(x,k-size(x,1)), 1+mod(k,2))]; endend

我對這個演算法比較心水，一方面容易理解，之前矢量化演算法從拉鏈矩陣對角線元素提取一步跳到sort的索引，加上裡面那個異常抽象的試湊構造矩陣，目瞪口呆之餘，氣氛很難熱鬧。相比之下，用diag逐條對角線處理更容易理解，另外一個好處是：無需考慮是否存在非0元素或全零對角線，甚至適用於不是方陣的情況。此外，還有一個小細節，「n<0」的判定被改成了「sign(n)」。

這個結構形式可改成arrayfun的形式，例如J.S.Kowontan的leading solution：

function ans = stop_zig(x,n) cell2mat(arrayfun(@(d)rot90(diag(rot90(x,sign(n)),d-length(x)),2*mod(d,2)),1:abs(n),uni,0));end

Pro44101. Adaptive ZigZag

前面兩個拉鏈矩陣的例子都比較形象，下面提到的這個問題，特徵則沒那麼明顯。

【鏈接】：https://ww2.mathworks.cn/matlabcentral/cody/problems/44101

【原題】：Unfold a 2-D matrix to a 1-D array in Adaptive zig-zag order, e.g., for matrix

[ 1 2 5 6; 3 4 7 8; 9 10 13 14; 11 12 15 16]

the resulting 1-D array should be

[ 1 5 9 13 2 3 4 6 7 8 10 11 12 14 15 16]

【解釋】：視輸入矩陣為4個分塊2×2子矩陣(可假定行列維度都為偶數)，依次提取每個子矩陣第1元素，匯總放在數組開頭，每個子矩陣剩下元素按第3、2、4，即主反對角線和最後一個元素的順序提取並順序擺放。

Test

Code Input

1

%%x =[1 2 5 6; 3 4 7 8; 9 10 13 14; 11 12 15 16] y_correct =[1 5 9 13 2 3 4 6 7 8 10 11 12 14 15 16];assert(isequal(your_fcn_name(x),y_correct))

2

%%x =1;y_correct = 1;assert(isequal(your_fcn_name(x),y_correct))

3

%%x =[1 2 5 6 1 6; 3 4 7 8 2 3; 9 10 13 14 6 8; 11 12 15 16 5 7] y_correct =[1 5 1 9 13 6 2 3 4 6 7 8 6 2 3 10 11 12 14 15 16 8 5 7];assert(isequal(your_fcn_name(x),y_correct))

【分析】：比較方便的做法是先決定怎樣提取分塊子矩陣，然後再各自內部提取相應元素。矩陣分塊提取可以用mat2cell函數，提取輸出結果放在cell數組內，例如：

>> a=randi(10,4,6)a = 9 7 10 10 5 7 10 1 10 5 10 1 2 3 2 9 8 9 10 6 10 2 10 10>> mat2cell(a,[2 2],[2 2 2])ans = 2×3 cell array {2×2 double} {2×2 double} {2×2 double} {2×2 double} {2×2 double} {2×2 double}>> ans{:}ans = 9 7 10 1ans = 2 3 10 6ans = 10 10 10 5ans = 2 9 10 2ans = 5 7 10 1ans = 8 9 10 10

mat2cell第2和第3參數分別是對原矩陣的分塊策略，顯然這段示例代碼和問題中要求一樣：2×2。

mat2cell、reshape和cellfun的命令組合

當看到mat2cell提取分塊子矩陣、解釋中提到的「分別處理」等等，應該能想到mat2cell後面順勢接cellfun用於處理子矩陣中元素的提取。

function ans = your_fcn_name(x) ans=x;[m,n]=size(x); if m>=2 reshape(mat2cell(x,2*ones(1,m/2),2*ones(1,n/2)),1,』』); cell2mat(cellfun(@(x)nonzeros(x),ans,un,0)); [ans(1,:) nonzeros(ans(2:end,:))]; endend

示例都是非0，因此用nonzeros簡化代碼，如果存在0元素，需要x(:)這類命令分行寫，這在下面的優化代碼中會繼續說明和介紹。同樣思路，彬彬的寫法略有不同：

function x = your_fcn_name(x)try t= cellfun(@(x)x([1,3,2,4]),mat2cell(x,repelem(2,size(x,1)/2),repelem(2,size(x,2)/2)),uni,0); t=cell2mat(t(:)); x=[t(:,1), reshape(t(:,2:end),1,)];end

一方面用repelem擴展全部維度為2的數組，並且結合輸入輸出共享變數x，用try流程保證可能出現的行列數小於2時不會錯誤輸出；此外注意try控制流程內第2行逗號表達式操作：

t=cell2mat(t(:));

自動把矩陣形式cell數組化成k×1形式傳給cell2mat處理，且每個cell中存儲的均是1×n數組，cell2mat自動將其排列為k×n矩陣。逗號表達式在MATLAB編程中應用非常豐富，有興趣可看看我和彬彬的書，有幾處結合實例對逗號表達式應用做了更詳盡的解釋。

J.S.Kowontan改變了cellfun的調用匿名函數句柄，進一步簡化上述代碼思路：

function x = your_fcn_name(x) try 0.5*size(x); cell2mat(cellfun(@(m)m(:),mat2cell(x,repelem(2,ans(2)),repelem(2,ans(1))),uni,0)); x = ans(union(1:4:numel(x),1:numel(x),stable));end

這段代碼有個隱藏在內的區別，以Test 1為例：mat2cell得到的是2×2cell數組，每個cell中由於冒號操作是4×1數組，沒有像之前彬彬用「t=cell2mat(t(:))」改變cell數組的排布方式，因而輸出4×4的double數組，而是在一個語句中一口氣套上cell2mat，形成的是8×2double數組：

>> x =[1 2 5 6; 3 4 7 8; 9 10 13 14; 11 12 15 16]x = 1 2 5 6 3 4 7 8 9 10 13 14 11 12 15 16

截止語句「cell2mat+cellfun(@(m)m(:),...)」的輸出：

ans=1 92 103 114 125 136 147 158 16

這個變化繼而影響到後續索引控制：用union命令提取，這同樣不容易想到，也異常精妙的索引操作。

>> t1=union(1:4:numel(x),1:numel(x),stable)t1 = 1 5 9 13 2 3 4 6 7 8 10 11 12 14 15 16

刻意構造出的8×2數組，需要按矩陣低維索引，每隔4個提取元素放在輸出數組最前面，其餘順次在後，使用union的stable參數，恰好滿足這個要求，這是個堪稱巧妙的構思，J.S.Kowontan對union函數似乎情有獨鍾，用union命令做構造不止一次在他的代碼中看到過。

玩魔方的permute+reshape的操控方式

高維數組維度變換中permute和reshape兩個函數應用普遍，不妨通過代碼熟悉用法。

function ans = your_fcn_name(x)b=numel(x);1;if b>1m=permute(reshape(permute(reshape(x,2,2,[]),[1,3,2]),2,2,[]),[2,1,3]);m=m(:);[m(1:4:b),m(mod(1:b,4)~=1)];endend

代碼如果出現permute，往往意味著3維，及3維以上的數組粉墨登場，這段代碼甚至連續應用2次permute+reshape的組合，不容易想到意圖，以Test1為例，拆開從裡到外逐層分析：

(1). 最裡層reshape(x,2,2,[])

注意到reshape函數除基本數組輸入x外，還有3個維度參數，說明4×4的矩陣x已正式通過reshape維度變換，成為2×2×4的3維數組；

>> t1=reshape(x,2,2,[])t1(:,:,1) = 1 9 3 11t1(:,:,2) = 2 10 4 12t1(:,:,3) = 5 13 7 15t1(:,:,4) = 6 14 8 16

(2). 觀察第1步所得3維數組中間變數t1，發現不滿足題意要求，所以用permute重新組織維度：

>> t2=permute(t1,[1 3 2])t2(:,:,1) = 1 2 5 6 3 4 7 8t2(:,:,2) = 9 10 13 14 11 12 15 16

簡單說明如何從t1變化到上述t2的：想像t1是一本4頁的書(維度3)，每頁都是2×2矩陣排布，現在換個方向看這個立的數組，從原來自頂而下視角換到側面自左向右看。permute函數第2參數[1 3 2]意味著列(第2)維度和頁(第3)維度互換，所以從原來的4頁變成2頁，從t1左側看過去，顯然是打亂t1，其每層第1列重新組合變成新數組t2的第1層，t1每層第2列變成t2的第2層。沒用過或者不熟悉permute函數的不妨參看幫助文件，此外我和彬彬書中對permute命令用法也結合實例做了詳細介紹。；

(3). 重新用reshape把數據換成2×2的多層形式：

>> t3=reshape(t2,2,2,[])t3(:,:,1) = 1 2 3 4t3(:,:,2) = 5 6 7 8t3(:,:,3) = 9 10 11 12t3(:,:,4) = 13 14 15 16

新變數t3又回到和第1步一樣的維度形式2×2×4，但元素位置有所變化，成功以「層」形式分塊。

(4). 最後用permute把3維數組t3每層元素轉置，便於下面步驟順序提取元素。

>> t4=permute(t3,[2 1 3])t4(:,:,1) = 1 3 2 4t4(:,:,2) = 5 7 6 8t4(:,:,3) = 9 11 10 12t4(:,:,4) = 13 15 14 16

這時發現只要把每層第1個元素提取出來，剩餘元素只需按低維索引順序導出即可，比如：

>> t4(:)ans = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

剩下的就比較容易了。J.S.Kowontan給出用union的改進版本

function ans = your_fcn_name(ans) try permute(reshape(permute(reshape(ans,size(ans,2),2,),213-0),2,2,),213-0); ans(union(1:4:numel(ans),1:numel(ans),stable));end

其中的』213』-『0』是縮減size的小手段，可以不去管它，直接用[2 1 3]即可。

這是個非常漂亮的多維數組從低維到高維再到低維的經典例子，也是不可多得學習reshape和permute在高維數組中應用的例子，多維數組元素在代碼編寫人手裡如同魔方，看似漫不經心、並且貌似重複的應用permute和reshape，只到元素排列完畢，再回溯倒推，發現來時的路幾乎步步殺機，而實際上卻履險如夷地完成了。

返璞歸真的基本命令組合

最後是Yurenchu完全用最基本命令組合寫出的代碼，所謂重劍無鋒、大巧不工，沒有花哨的函數，也沒有叫人困惑的演算法，從頭到尾和索引的構造在打交道，明面上唯一用到的函數只有人人熟悉的size，Yurenchu的MATLAB代碼基本功實在令人嘆服。如果仔細查看if流程中的第2和第4句，你會發現，令人嘆服的MATLAB代碼基本功，實際上是建立在其對矩陣元素之間內在聯繫的敏銳嗅覺上的。

function x = your_fcn_name(x) [a,b] = size(x); if a>1 && b>1 b = (1:2:a) + (1:2:b)*a - a; b = b(:); a = [a; 1; a+1]+b; x = x([b a(:)]) endend

小結

本篇主要討論在矩陣中如何尋找和按照特定順序提取具有一定規律的元素，這類問題常見於Cody，我和彬彬書里也討論了類似bulleye matrix這樣的問題，本篇所討論內容雖然從表面形式和要求上，迥異於書中的bulleye matrix，本質上仍然考驗代碼編寫人對矩陣元素內在規律特徵的把握，以及尋找合適命令組合實現的能力。因此，多維數組命令也好，基本矩陣索引操控也罷，或具有特殊形式的矩陣構造，都是這種本質能力在軟體應用時的實際反映。

上述問題代碼的命令可用下圖表示：