CodyNote014:再談動態長度子序列(Part.3)

08-25

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

馬良，祁彬彬

序言

動態長度子序列問題對編寫MATLAB代碼培養開闊思路，深刻了解各類控制流程和命令組合有很大好處，同時又具有相當難度，需全力以赴絞盡腦汁，當提交答案在世界級舞台上供人閱覽，偶爾獲得認同時，成就感自不必提！

下面將繼續介紹幾個動態子序列的問題，慢慢會發現，這類問題代碼的高度靈活與開放性，充分激發了Cody幾位超級高手的征服欲，往往竭盡所能，提供優秀經典的1種甚至多種代碼方案，不寫出令自己滿意的最優化方案根本不罷休。鑽研這類問題，是全面掌握MATLAB編程技巧的最優途徑、是遭遇MATLAB代碼編寫的瓶頸後二度突破的通天之路。而且屬於自己閉門造車，悶頭苦學任何教材無法獲得的寶貴提高機會。例如CodyNote013中，Alfonso給出的最小二乘解方案，從MATLAB視角說是演算法、數學建模、矩陣理論和MATLAB編程語言四位一體無縫結合的經典佳作，也毫不過分。看完這種代碼，再去學習普通教材里的標準範式代碼流程，就往往能有不一樣的理解層次，而我認為，創新的火花，就隱藏在這種無法言喻的莫名契機中。

我在和彬彬的書中提到這樣一句話：「寫出優秀代碼前，首先你得能夠欣賞優秀代碼」。

創新思維和能力培養，蓋不過如是。

特徵子序列判斷與提取

下面的數據最大跨度求解問題，求解時感覺印象異常深刻，因為最開始我用mode走了一大截彎路。

Pro317. Find the stride of the longest skip sequence

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

【原題】：We define a skip sequence as a regularly-spaced list of integers such as might be generated by MATLABs colon operator. We will call the inter-element increment the stride. So the vector 2:3:17 or [2 5 8 11 14 17] is a six-element skip sequence with stride 3.

Given the vector a, your job is to find the stride associated with the longest skip sequence you can assemble using any of the elements of a in any order. You can assume that stride is positive and unique.

Example:

input

a = [1 5 3 11 7 2 4 9]

output stride is 2

since from the elements of a we can build the six-element sequence [1 3 5 7 9 11].

【解釋】：出題人Ned Gulley似為Mathworks公司專門負責Cody運轉維護的工程師，此人也是出Cody常規96題的Cody Team成員之一。這題對求解人演算法水平有著不低的要求，題目立意不落窠臼，「教材流」的兄弟大概率被折磨得死去活來。大意是給定大小順序打亂的行數組，從中找出長度最長等差序列的公差項。例如題目示例中的input，其中所有元素能組成的等差序列中，最長是1:2:11，所以公差值為2。

Test

Code Input

%%a = [1 5 3 11 7 2 4 9];stride = 2;assert(isequal(skip_sequence_stride(a),stride))

%%a = [1:5:20 23:3:42 2:9:100];stride = 9;assert(isequal(skip_sequence_stride(a),stride))

%%a = [2:2:22 13:17];a = a(randperm(length(a)));stride = 2;assert(isequal(skip_sequence_stride(a),stride))

%%a = 37:5:120;a = a(randperm(length(a)));stride = 5;assert(isequal(skip_sequence_stride(a),stride))

%%a = [1:5 101:10:171 201:205];a = a(randperm(length(a)));stride = 10;assert(isequal(skip_sequence_stride(a),stride))

%%a = [7:17:302 primes(300)];a = sort(a);stride = 17;assert(isequal(skip_sequence_stride(a),stride))

【分析】：問題看似無從下手，常規控制流程例如for、if等等，都不知道應該從哪個位置開始寫起。如果歸攏思路，會發現難點無非在於：怎麼能枚舉該輸入行數組所能構造出的各個等差序列。

一個錯誤思路引發的慘案

這事兒說起來比較噁心，我開始直覺可以用隱式擴展自相減，因為找最長等差序列公差項，分別相減後，公差項似乎就是相減後矩陣里，出現頻次最高的數。當然，默認公差項不為零，Test Suite也證實了這一點，於是興沖沖按這個思路擴展，得到的代碼非常簡單，就兩句：

a=unique(a);m=nonzeros(triu(a-a));

不管對錯，先看代碼：

(1). 第1句中unique排除重複項轉換為升序排列，這裡要注意：unique、union這些命令都具備隱含排序能力，能對輸出結果自動排序，以test1為例：

>> a = [1 5 3 11 7 2 4 9];>> unique(a)ans = 1 2 3 4 5 7 9 11

容易看出，unique不帶任何後置參數，排除重複元素後的結果自動升序，很多場合下用得較普遍。當然有的情況下也需要元素仍按原順序排布，那就用到後置的「setorder」參數，這個參數共2個選項：sorted(default)和stable，選擇後一個參數「stable」就按原序列順序排布：

>> b=[1 3 12 3 7 4 7 1 212 33 12 33]b = 1 3 12 3 7 4 7 1 212 33 12 33>> unique(b,stable)ans = 1 3 12 7 4 212 33

unique命令後置參數設置方法簡單，但很靈活，比如還有按整行是否相同進行剔除的rows參數，按出現次序排除重複的first和last，即：到底剔除先出現還是後出現的元素，詳見help。

重點是第2句，先說用到的mode+nonzeros組合，mode在NoteCody013中已經介紹：找到最大頻次元素，nonzeros是2006a版本以後的稀疏矩陣命令，用於檢索矩陣非零元素，可以把矩陣全部非零元素按列形式排布提取出來，我習慣在求和時替代sum(sum(...))的寫法。

再說裡層的「triu(a-a』)」，注意：此時a是經unique剔重的升序列。雖探討多次，但還要再重申單一維度擴展，也就是隱式擴展的重要性，如果用R2016b後的版本，應對此用法引起足夠重視，即使不會用，也要理解它可能得到的結果：行序列a減其自身轉置，這有悖於普通線性代數教材中的提法，但相當多情況下，卻屬於便利的操作方式。數組(1×n)-數組(n×1)，儘管不同維，但減數列維度為1，相當於被減數分別減去後面n×1數組每個元素，排列成n行，結果是n×n方陣。

>> a=unique(a)a = 1 2 3 4 5 7 9 11>> a-aans = 0 1 2 3 4 6 8 10 -1 0 1 2 3 5 7 9 -2 -1 0 1 2 4 6 8 -3 -2 -1 0 1 3 5 7 -4 -3 -2 -1 0 2 4 6 -6 -5 -4 -3 -2 0 2 4 -8 -7 -6 -5 -4 -2 0 2 -10 -9 -8 -7 -6 -4 -2 0

按前面分析從答案逆推：既然1:2:11為最長等差序列，2這個元素相減過程中似乎，注意，是「似乎」出現頻次最高就應該才對，下三角陣相當於上三角陣相反數，triu命令把下方元素置零，nonzeros提取全部非零元素，用mode得最大頻次數，即：公差項，看結果對不對：

>> mode(nonzeros(triu(a-a)))ans = 2

這個結果是對的！繼續測試第2、3和第4個test的assertion的結果，也都吻合，這讓我很高興，但接著代入第5和第6個test開始出現錯誤，測試5結果應當為10，而運行結果為1：

>> a = [1:5 101:10:171 201:205];>> a=unique(a);mode(nonzeros(triu(a-a)))ans = 1

測試6正確結果應為17，運行上述代碼結果為30：

>> a = [7:17:302 primes(300)];>> a=unique(a);mode(nonzeros(triu(a-a)))ans = 30

因此mode+triu+單一維擴展相減的思路行不通，至少在某些地方做出修正前是行不通的。直接用test5來說明前述代碼：「mode(triu(unique(a)-unique(a)』))」的錯誤，以test5序列為例：

a = [1:5 101:10:171 201:205];>> aa = 1 2 3 4 5 101 111 121 131 141 151 161 171 201 202 203 204 205

顯然，公差為1的子序列被分成前後兩部分：開始的1:5和最後的201:205，儘管兩部分序列總長度最大( $5+5=10$ )，但兩段不能合併組成完整公差項為1的連續序列，中間101:171儘管總長度 $8<5+5=10$ ，但分別大於2個公差項為1的子序列，的確是最長連續子序列。

因此，test5正確結果不是1，應該是10。前述代碼雖然短，但未考慮到分段的情況用mode盲目合併，得到錯誤結論。而且堅持不放棄mode+triu+單一維擴展，個人感覺轉回來非常困難，當時費了好大勁兒，花了很長時間，硬掰出個很噁心的代碼，後續又是histcounts，又是ismember，還有regexp正則查找，加上cellfun嵌套arrayfun，兵兵乓乓弄得好不熱鬧！可事實卻是代碼的效率低下且冗長乏味，在低級趣味、裝B失敗的道路上漸行漸遠：

function ans = skip_sequence_stride(a)a=unique(a);m=nonzeros(triu(a-a));[t,bin]=histcounts(m,[-inf;unique(m)]);temp=[bin(2:end-1);t(2:end)];stridek=temp(1,temp(2,:)>=median(temp(2,:)));arrayfun(@(n)max(arrayfun(@(k)max(cellfun(@numel,regexp(num2str(ismember(a(k):n:a(end),a),-6),1+,match))),1:numel(a))),stridek);stridek(ans==max(ans));end

儘管問題勉強通過全部test suite，不得不承認，這屬於典型的失敗代碼案例。這是思路發生偏差，照死里猛鑽牛角尖導致的慘案，也許通過一些技巧能把代碼內部再做一些小修小補的優化，但個人認為不值得，因為有更好的方案。

一般的循環判斷

循環判斷，其實也比前面這個「偽」矢量化代碼要好得多：

function ans = skip_sequence_stride(a)n=0;for m=a for u=a for i=1:20 k=m:i:u; if ismember(k,a)& length(k)>n n=length(k);ans=i; end end endend

思路灰強樸實：三重循環在a數組內，以外層m為起點，內層u為末端點逐步前移，試探構造子序列「m:i:u」，尋找不同長度組合下的連續等差數列，用ismember判定數列「m:i:u」中元素是否全部為a中元素，並確定該子序列是否長度大於前一子序列長度n，當然，n的長度初值為0，找到合適子序列後要在裡層if判定更新一下。

利用圖論Dulmage-Mendelsohn 分解命令dmperm

Alfonso用圖論工具箱命令dmperm，把基本3重循環降為一重循環，簡單搜索發現有時dmperm函數應用在二值圖像連通域判斷，如bwlabel函數使用Dulmage-Mendelsohn 分解消除等價對，鏈接(https://blog.csdn.net/muye_CSDN/article/details/77447958?locationNum=7&fps=1 )給出以深度優先遍歷替換等價對的等效做法，該函數全格式調用方法如下：

[p,q,r,s,cc,rr] = dmperm(A)

(1). 如果列 j 與行 i 匹配，p = dmperm(A) 得到結果為向量 p，p(j) = i：

如列j與其不匹配，結果為0。如A是具有完整結構秩的方陣，則p是最大匹配行置換且A(p,:)包含非零對角線。A 的結構秩是sprank(A) = sum(p>0)。

(2). [p,q,r,s,cc,rr] = dmperm(A)（其中A無需是方陣或完整結構秩）計算 A 的 Dulmage-Mendelsohn 分解。p 和 q 分別是行和列置換向量，A(p,q) 包含分塊上三角形。r 和 s 是索引向量，指示精細分解的塊邊界。cc 和 rr 是長度為 5 的向量，指示粗略分解的塊邊界，所謂粗略分解(coase decomposition)指C=A(p,q)被粗略分成4×4粗略塊，因此行、列邊界(cc,rr)為4+1=5，例如：

>> x=randi([0 1],4)x = 0 0 1 0 1 0 1 1 0 0 0 0 0 1 0 0>> [p,q,r,s,cc,rr]=dmperm(x)p = 2 4 1 3q = 4 1 2 3r = 1 2 3 4 5s = 1 3 4 5 5cc = 1 2 3 5 5rr = 1 2 4 4 5

再搜wiki百科(https://en.wikipedia.org/wiki/Dulmage%E2%80%93Mendelsohn_decomposition )，找到Dulmage-Mendelsohn分解是圖論中，用於分區二分圖(bipartite graph)所有頂點至各個子集合的最大匹配。相鄰頂點在同一子集。如果不學習圖論，可以略過不看，直接跳到Alfonso提交的代碼：

function ans = skip_sequence_stride(a) d=a-a; for n=unique(d) [~,~,r]=dmperm(~d|abs(d)==n); ans(max(diff(r)))=n; end ans(end);end

這段代碼：對邏輯矩陣「~d|abs(d)==n」，按Dulmage-Mendelsohn 分解，所得分解塊矩陣C的精細分塊行索引向量r的數值。第1行d=a-a』為單一維擴展，詳見第1小節我的錯誤解答，但那時用得不對，Alfonso在這裡用dmperm命令使其重煥新生。以test1為例，循環體對d剔重得到結果為-10:10：

>> n=unique(d)n = -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10

進入循環後，當n<=0，即-10:0這階段，dmperm輸入邏輯矩陣「~d|abs(d)==n」的後半段因n<=0，顯然是全0陣，「~d|abs(d)==n」此時結果實際上是「~d」，也就是單位邏輯矩陣，因此11次循環結果的r是1:9。當循環進入n>=1後，後半段不再是全0陣，Dulmage-Mendelsohn 分解結果開始變化，頂點子集精細分塊的行索引r不再是1:9，例如當n=2時，有t=~d|abs(d)==n，結果為：

t = 8×8 logical 數組 1 0 1 0 0 0 0 0 0 1 1 0 1 0 0 0 1 1 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 1 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 1 1 0 0 0 0 1 1 0 0 1

精細分塊行索引r的結果相應為：

>> [~,~,r]=dmperm(~d|abs(d)==n)r = 1 7 9

按照循環體最後一行「ans(max(diff(r)))=n」，此時ans第6項應為n，也就是2，其餘類推，注意這一行中的max(diff(r))，也就是索引跨度越大，在ans中的位置就應當越靠後，循環結束，取最後一項就是最長等差數列的跨度值。

ismember命令的正確打開方式

在第1小節我用ismember+regexp檢索最長等差序列，但用起來磕磕巴巴相當糟糕，那個味道和感覺全然不對。下面看ismember命令在本問題中的正確打開方式。

function stride = skip_sequence_stride(a) [b,c]=ndgrid(sort(a),1:max(a)); [~,stride]=max(max(arrayfun(@(x,y)find(~ismember(x:y:max(a)+y,a),1),b,c)));end

相比前面的圖論命令dmperm，這個演算法好理解得多，它其實是之前三重循環演算法的MATLAB矢量化壓縮版。從第1句的網格點布點函數ndgrid輸入參數其實就能看出端倪：這是要把所有段落的等差序列用「包絡」的形式枚舉出來，話雖然講起來非常輕鬆，可是真要實現，尤其是被壓縮到1-2行中實現出來，卻絕對不是什麼輕鬆的事情，尤其是應用arrayfun這一句中，包含非常有意思，在理解上容易形成誤區的arrayfun用法。仍以test1中的變數a為例：

>> aa = 1 5 3 11 7 2 4 9

用代碼第1句ndgrid對其進行網格布點，得到b和c都是8×11的矩陣，此時注意第2句中，ndgrid輸出b和c直接放在arrayfun中，所說誤區是指本來按arrayfun字面意思，顧名思義應只與數組array有關，所以常在arrayfun第2輸入參數中使用諸如：1:numel(x)、1:length(x)做循環節，矩陣做輸入也不是不可以，可往往都是像普通循環那樣，整列讀取，但這個例子就讓我開了眼界，不妨去掉內部第1參數中的find、ismember，直接看看「x:y:max(a)+y」能夠在arrayfun以b和c做直接調用的情況下，形成什麼輸出結果：

>> ty=arrayfun(@(x,y)x:y:max(a)+y,b,c,un,0)ty = 8×11 cell 數組 1 至 10 列 {1×12 double} {1×7 double} {1×5 double} {1×4 double} {1×4 double} {1×3 double} {1×3 double} {1×3 double} {1×3 double} {1×3 double} {1×11 double} {1×6 double} {1×5 double} {1×4 double} {1×3 double} {1×3 double} {1×3 double} {1×3 double} {1×3 double} {1×2 double} {1×10 double} {1×6 double} {1×4 double} {1×4 double} {1×3 double} {1×3 double} {1×3 double} {1×3 double} {1×2 double} {1×2 double} {1×9 double} {1×5 double} {1×4 double} {1×3 double} {1×3 double} {1×3 double} {1×3 double} {1×2 double} {1×2 double} {1×2 double} {1×8 double} {1×5 double} {1×4 double} {1×3 double} {1×3 double} {1×3 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×6 double} {1×4 double} {1×3 double} {1×3 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×4 double} {1×3 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} 11 列 {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double} {1×2 double}

結果ty是和原來b和c同維的cell數組，以ty第1個元素為例，這是個1×12的double數組，它是怎麼來的呢？

>> ty{1,1}ans = 1 2 3 4 5 6 7 8 9 10 11 12>> b(1):c(1):max(a)+c(1)ans = 1 2 3 4 5 6 7 8 9 10 11 12

是的，並沒有像我之前以為的那樣，arrayfun就只能處理1×n或者n×1的數組，它也是可以點對點操作的，所以其他問題就比較好理解了：在矩陣b以內，循環前移，錯位以c(m,n)為步距獲得等差序列，再右移c(m,n)，用find+~ismember判斷第1個不屬於a序列的元素位置，例如ty{1,1}，按照find(~ismember(b(1):c(1):max(a)+c(1)),1)得到的第1個非a元素位置是6，說明前5個都是，其餘類推。因此在最終獲得的8×11矩陣中，囊括的是所有不同位置移位、構造等差序列中第1個不屬於a的元素位置，以test1為例：

>> t1=arrayfun(@(x,y)find(~ismember(x:y:max(a)+y,a),1),b,c)t1 = 6 7 4 4 2 3 2 3 2 3 2 5 3 3 2 3 2 3 2 3 2 2 4 6 2 4 2 3 2 3 2 2 2 3 2 3 2 3 2 3 2 2 2 2 2 5 2 3 2 3 2 2 2 2 2 2 4 2 3 2 2 2 2 2 2 2 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

外層第1個max獲得每列中的這個非a元素的最大首索引的數值：

>> t2=max(t1)t2 = 6 7 4 4 3 3 3 3 3 3 2

最外層max不是要最大值，而是該最值的索引位：顯然t2最大元素是7，所在位置2，這是最終計算結果。

這個結果看起來好像略懸乎，為放心起見，逆推結果產生的過程，就能看出結果2的出處。

第1階段逆推，由矩陣t1，全矩陣最大值為7，發生在第(1,2)索引位，因此把(1,2)代入得：

>> b(1,2):c(1,2):max(a)+c(1,2)ans = 1 3 5 7 9 11 13

前6個元素都是a中元素，只有第7個13不是，所以1:2:11就是在a中能構造出的最長等差序列。

>> ~ismember(ans,a)ans = 1×7 logical 數組 0 0 0 0 0 0 1>> find(ans)ans = 7

第2階段逆推，判斷索引(1,2)怎麼來的？回到之前的ndgrid：

>> [b,c]=ndgrid(sort(a),1:max(a))b = 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 7 7 7 7 7 7 7 7 7 7 7 9 9 9 9 9 9 9 9 9 9 9 11 11 11 11 11 11 11 11 11 11 11c = 1 2 3 4 5 6 7 8 9 10 11 1 2 3 4 5 6 7 8 9 10 11 1 2 3 4 5 6 7 8 9 10 11 1 2 3 4 5 6 7 8 9 10 11 1 2 3 4 5 6 7 8 9 10 11 1 2 3 4 5 6 7 8 9 10 11 1 2 3 4 5 6 7 8 9 10 11 1 2 3 4 5 6 7 8 9 10 11

原來網格布點包絡是保證完全不遺漏提取所有與a元素有關序列的關鍵，在b(1,2):c(1,2):max(a)+c(1,2)中最後加的c(1,2)也是為了保證尾部元素不遺漏。

Alfonso這個解法和思路真是不要太拽！

當然，使用arrayfun還有其他同樣出色的方案，例如@bmtran增加一個arrayfun但省去ndgrid：

function stride = skip_sequence_stride(a)[l,stride] = max( arrayfun( @(s)max( arrayfun(@(e)find(~ismember(e:s:max(a)+s,a),1), a) ), 1:max(a) ) );end

Alfonso在原來基礎上對代碼又做了兩種優化方案，如下是第1種：

function stride = skip_sequence_stride(a)[l,stride] = max(arrayfun(@(s)norm(sign(expm(+bsxfun(@eq,a,a+s))),1),1:max(a)));end

這裡面同樣使用了和圖論理論有關的命令組合：「norm(sign(expm))」，在另外一個關於求解有向圖中連通域個數的問題(https://ww2.mathworks.cn/matlabcentral/cody/problems/2372 )中，Alfonso同樣用到了類似的組合，比如下面的這個矩陣x：

x=[0 1 0 0 0 0 ; 1 0 0 0 0 0; 0 0 0 1 0 0; 0 0 1 0 0 0 ; 0 0 0 0 0 1; 0 0 0 0 1 0];

如果用digraph命令生成有向圖G，並plot得到如下結果：

>> x=[0 1 0 0 0 0;1 0 0 0 0 0; 0 0 0 1 0 0;0 0 1 0 0 0 ;0 0 0 0 0 1;0 0 0 0 1 0];>> G=digraph(x);>> plot(G)

圖1 矩陣x的有向圖

顯然，封閉連通域個數應為3，可通過MATLAB2015b中新增的圖論工具箱命令conncomp求解：

>> max(conncomp(digraph(x)))ans = 3

Alfonso提交答案沒用圖論工具箱命令，因為提交時間為2014年，這個工具箱還沒出現，他的做法是：

>> rank(sign(expm(x)))ans = 3

這從側面體現了大神Alfonso在矩陣理論方面的洞察力和對MATLAB函數極其全面和深刻的理解。同時也提醒我們sign+expm函數組合與圖像連通域判定有著比較密切的理論聯繫。

上述解答並非全劇終， Alfonso接著又祭出矩陣多項式求值命令polyvalm+norm的組合大招，得到該解法另一種變體，同時代碼size相比之前更小：

function stride = skip_sequence_stride(a)[l,stride] = max(arrayfun(@(s)norm(polyvalm(a,a==a+s),1),1:max(a)))end

應該都知道多項式和矩陣特徵值間存在密切聯繫，簡單舉個例子：

>> A = [1 2 3; 4 5 6; 7 8 0]A = 1 2 3 4 5 6 7 8 0>> poly(A)ans = 1.0000 -6.0000 -72.0000 -27.0000>> roots(ans)ans = 12.1229 -5.7345 -0.3884>> eig(A)ans = 12.1229 -0.3884 -5.7345

根據線性代數知識，poly得到矩陣特徵多項式，其根恰為該矩陣的特徵值。儘管都是多項式求值命令，儘管調用方式類似，但polyvalm命令與polyval的不同之處主要是前者遵從矩陣整體運算方式，區別如下：

>> p=1:3;A=eye(2);>> t1=polyval(p,A)t1 = 6 3 3 6>> A.^2+2*A+3ans = 6 3 3 6>> t2=polyvalm(p,A)t2 = 6 0 0 6>> A^2+2*A+3*eye(2)ans = 6 0 0 6

問題在於矩陣多項式求值命令描述矩陣基本特徵的方式。上述運算結果可知，它與polyval存在差異：後者是單值的項式求值代入，而前者則需要引入範數概念。查閱矩陣論教材，發現範數基本概念其實是一般線性空間中，刻畫收斂性質的廣義長度。這個長度數值的概念，與問題中的「最長」等差序列有關，該長度值和norm，即範數數據間，存在正相關的聯繫，正是基於這個概念的直覺，Alfonso提出用norm評價polyvalm對向量a，在「a==a+s」矩陣處的特徵多項式求值結果的概念。為方便理解，我把Alfonso用arrayfun簡化的命令，重新改寫成for循環格式：

function ans=skip_sequence_stride(a)for i=1:max(a) ans(i)=norm(polyvalm(a,a==a+i))end[~,ans]=max(ans);end

以test1為例，得到的結果如下：

>> skip_sequence_stride(a)ans = 21.6746ans = 21.6746 25.1589ans = 21.6746 25.1589 12.7009ans = 21.6746 25.1589 12.7009 12.7009ans = 21.6746 25.1589 12.7009 12.7009 11.2195ans = 21.6746 25.1589 12.7009 12.7009 11.2195 11.2195ans = 21.6746 25.1589 12.7009 12.7009 11.2195 11.2195 11.2195ans = 21.6746 25.1589 12.7009 12.7009 11.2195 11.2195 11.2195 11.2195ans = 21.6746 25.1589 12.7009 12.7009 11.2195 11.2195 11.2195 11.2195 11.2195ans = 21.6746 25.1589 12.7009 12.7009 11.2195 11.2195 11.2195 11.2195 11.2195 11.2195ans = 21.6746 25.1589 12.7009 12.7009 11.2195 11.2195 11.2195 11.2195 11.2195 11.2195 9.0000ans = 2

Aofonso這個解法最令人印象深刻的還不是他想到polyvalm+norm組合能夠等效取代norm(sign(expm))——當然，已經夠嚇人了。我意思是，在polyvalm內部的矩陣構造「polyvalm(a,a==a+s)」，以及用norm求取其1範數的結構設計，才是讓我頭頂冒油汗的真正原因：這要有多麼強悍的矩陣思維，才能在如此短時間內，想到這麼精巧的思路，並將其一一解構、替換為合理的代碼？

完全的矢量化：shiftdim+cumprod+ismember+max

最後是另一位大神LY Cao的答案，他對MATLAB函數命令的理解，與Alfonso相比不遑多讓，此外，這是所有代碼中，矢量化最徹底的——連arrayfun都去掉了：

function k = skip_sequence_stride(a)[~,k] = max(max(sum(cumprod(ismember(a+find(a).*str2num(shiftdim(1:100,-1)),a)))));end

先說說命令shiftdim，用於數組內部元素的維度變換，也就是說，它可以把數據內的數據維度進行輪轉，例如一個簡單的例子：

>> arrayfun(@(i)shiftdim(1:3,i),-2:4,un,0)ans = 1×7 cell 數組 {1×1×1×3 double} {1×1×3 double} {1×3 double} {3×1 double} {1×3 double} {3×1 double} {1×3 double}K>> ans{:}ans(:,:,1,1) = 1ans(:,:,1,2) = 2ans(:,:,1,3) = 3ans(:,:,1) = 1ans(:,:,2) = 2ans(:,:,3) = 3ans = 1 2 3ans = 1 2 3ans = 1 2 3ans = 1 2 3ans = 1 2 3

把數組1:3，其維度從-2變化到4，一共7次變換，第1次按-2變換，結果變成4-D數據，這個數據的來源是，1:3其維度按size命令可知為1×3，也就是[1 3]，雖然後面3維和4維甚至更高維度都沒寫，默認為1，即：[1 3 1 1]；-2代表維度從末尾向前一共前移2位，也就是：「[1 3]→[1 1 3 1]→[1 1 1 3]」，當轉換維度參數變正值後，因最小維度只能為1，此時[1 3]已經沒有「0維度」可轉換，第3和第4次又轉回原來[1 3]和[3 1]維度。二維或者更高維度數據也是如此：

>> t=[1 2 3;4 5 6];K>> arrayfun(@(i)shiftdim(t,i),-2:4,un,0)ans = 1×7 cell 數組 {1×1×2×3 double} {1×2×3 double} {2×3 double} {3×2 double} {2×3 double} {3×2 double} {2×3 double}

因此，LY Cao代碼中的「shiftdim(1:100,-1)」其結果維度應該是3維數據1×1×100。緊接著，與之做點乘的find(a)』相當於(1:numel(a))』，以test1為例，這是8×1向量，和之前3位數據1×1×100點乘，根據維度對應操作的基本原理，得到8×8×100的3維數組，其實是把1:100每個數對應和(1:8)』相乘，並存儲在3維數組的第3維度上，講到這裡，為什麼能省略arrayfun命令就很清楚了，因為矢量化運算被放在第3維度。上述從MATLAB編程角度，即數據維度協調說明問題。如果從問題視角，這種操作目標就是枚舉和構造一堆不同的等差數列，用ismember(...,a)來分別比較這個複雜的8×8×100數組中，哪些是a中的元素，哪些不是。

下面又是亮點：為何使用cumprod命令？仍以test1中a(1×8)為例，因find(a)』對1:100每個數都乘了一遍，給出從1到100共100個公差項，且a每個數值都有可能從8×8矩陣第1行起頭，使用cumprod命令對ismember獲得的邏輯數組進行累積乘積，最長連續1就一定是最大等差數列長度，因為ismember返回結果只有0和1，累積乘積不為零就必須連續不是1，如果乘到0，那麼後續所有數字都應當是0.剩下的就是用sum把所有列上的1累加，最多的就是最長等差序列，連續兩個max和之前敘述相同，不再贅述。

從低維數組轉換為高維數組的過程連貫而一氣呵成，避開for循環，一次完成對數據的乘運算，同樣是難得的精品代碼。

當然，矢量化運算，高維數組非常吃內存，數據比較大的時候，效率可能反而不如三重循環，但從學習角度，毫無疑問這個代碼對高維數組運算以及高、低維數組轉換操作的理解具有極大好處。

小結

這篇內容只探討1個問題，信息量卻絲毫不遜色前一篇，內容中包含矩陣論、圖論、高維數組運算等內容的綜合，從題目求解多種代碼，實錘了對動態長度子序列問題的代碼編寫，需要對矩陣運算基礎理論、MATLAB基礎矩陣和多項式運算命令的理解深刻到一定高度。

沒錯！歡迎來到MATLAB編程的深水區！你看到的，就是實打實、世界最頂級MATLAB高手全力以赴，或者接近火力全開狀態下表現出的水準。