CodyNote002：富有彈性的數組長度(part.2)

08-21

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

馬良，祁彬彬

接前一部分：

Pro42591.Produce the following matrix

本節是另一類彈性長度矩陣的問題，要求構造2維矩陣，其行或列之間具有顯著變化規律特徵。

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

【原題】：Produce the following matrix

x = [2 3 4] y_correct = [1 1/2 1/3; 2 1 1/4; 3 4 1];

【解釋】：題目意思都在例子中，實際上就是矩陣以主對角線分界，對稱元素互為倒數，主對角線全是元素1。

Test

Code Input

%%filetext = fileread(matrixManipulation.m);assert(isempty(strfind(filetext, system)))

%%x = 1;y_correct = ones(2);assert(isequal(matrixManipulation(x),y_correct))

%%x = 2;y_correct = [1 1/2; 2 1];assert(isequal(matrixManipulation(x),y_correct))

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

%%x = 2:7;y_correct = [1 1/2 1/3 1/4; 2 1 1/5 1/6; 3 5 1 1/7; 4 6 7 1];assert(isequal(matrixManipulation(x),y_correct))

【分析】：題目意思看似直白而簡單，實則不太容易，常規思路要先獲得下三角陣元素，轉置求倒數，二者相加。所以不妨看看各路高手是怎麼處理這個問題的，其中包含了不少實用的編程技巧。

for循環直譯的方式

下面給出一段平鋪直敘的代碼作為分析和述評的開始。

function y = matrixManipulation(x)count=1;n=round((sqrt(1+8*length(x))+1)/2);y=eye(n);for i=1:n for j=i+1:n y(j,i) = x(count); y(i,j) = 1/x(count); count=count+1; endendend

代碼難點在第1句：矩陣維度計算給出了計算公式：

n=round((sqrt(1+8*length(x))+1)/2)

用eye(n)得到全1主對角線後，以二重循環逐次填充。

利用tril+find完成索引查找

如果要求不是太嚴格，題目的test suite要求的2維矩陣的維度都不大，可以直接用cumsum枚舉：

>> cumsum(1:10)ans = 1 3 6 10 15 21 28 36 45 55

接著從ans中找到與題目輸入數組長度+1相等的數值所在索引位位置，這就是輸出矩陣的維度。

>> find(ans==numel(2:7))+1ans = 4

所以test7中，輸入x=2:7，輸出填充矩陣的行列數應為4.總體代碼：

function ans = matrixManipulation(x)find(cumsum(1:10)==numel(x))+1;t=tril(ones(ans),-1);[ik,jk]=find(t);for i=1:numel(x), t(ik(i),jk(i))=x(i);endtriu(1./(t+~t))+t;end

最後利用上三角和下三角元素間的倒數關係相加。

矩陣翻轉和邏輯索引

function y = matrixManipulation(x)str2num [2 2 3 3 3 4];y=ones(ans(length(x)));y(tril(1|ones(ans(length(x))),-1))=x;y=1./fliplr(rot90(y,-1));y(tril(1|ones(ans(length(x))),-1))=x;end

呃...第1行的維度枚舉方式未免太...「量身定做」了一點兒吧？直接就枚舉了...

先不管它，這段代碼的關鍵就是直接通過構造符合要求的邏輯索引進行賦值，相比前一種解法中通過循環賦值，在MATLAB中，這才是值得推薦的矢量化代碼技巧。此外，也要注意邏輯數組構造的方式採用了「1|ones(...)」的形式，而不是用logical命令，這是cody中縮減代碼size的常用技巧。

還有一種解法綜合了第1和第3種解法的優點：

function y = matrixManipulation(x) N = (1+sqrt(1+8*length(x)))/2; y = eye(N); mask = tril(true(N),-1); y(mask) = x; y = y; y(mask) = 1./x; y = y;end

去掉flip+rot90的翻轉以轉置代替，維度計算以通用方法，邏輯索引選擇使用true函數。

關於維度的計算，Tim提供了另一種方式：

function y=matrixManipulation(x)y=ones(max(2,0.5+sqrt(2*max(x)-1.75)));k=find(tril(y,-1));y(k)=x;y=y;y(k)=1./x;y=y;end

注意，矩陣維度公式「0.5+sqrt(2*max(x)-1.75)」只適合n>2的情況，所以用max保證小於2的計算結果能用2「兜」住而不至於出錯。這和前一道題目提到的max用法一致。維度的計算還可以用roots直接計算：

function y = matrixManipulation(x)n = max(roots([1,-1,-2 * length(x)]));y = triu(ones(n));y(~y) = x;y = y - 1 + y.1 ;end

或者ceil+sqrt的組合：

function y = matrixManipulation(x) n = ceil(sqrt(numel(x)*2)); y = eye(n); y(~triu(true(n))) = x; y = max(y,triu(1./y.));end

Jan Orwat這段代碼還有個隱藏的精彩構思：最後一行用max在y和倒數間做元素對位比較。我和彬彬的書中也有幾個題目介紹了max/min在題目的類似用法，這個做法省去了相加或者分別按邏輯索引賦值的運算。基於序言中提及的彈性擴維概念，eye的賦值也同樣能夠省略：

function y = matrixManipulation(x) y = ones(ceil(sqrt(numel(x)*2))); mask = ~triu(y); y(mask) = 1./x; y=y.; y(mask) = x;end

再看看劉鵬的代碼:

function y = matrixManipulation(x)ones(ceil(sqrt(2*numel(x))));ans(~triu(ans)) = x;y = ans./ans.;end

第1行通過ceil+sqrt確定矩陣維度，並構造全1陣，這和其他代碼類似;

第2行的構思則異常精巧，首先，其本身考慮到了triu函數自帶主對角線的特徵，通過取反邏輯操作獲得不含主對角線的下三角邏輯索引，由於下三角的元素數量和輸入x正好相同，因此按列排布方向把x的元素逐一賦值給矩陣的下三角部分，這是索引矢量化的經典代碼，非常值得學習：

>> ~triu(ones(5))ans = 5×5 logical array 0 0 0 0 0 1 0 0 0 0 1 1 0 0 0 1 1 1 0 0 1 1 1 1 0>> x=2:11;>> t=ones(5);>> t(ans)=xt = 1 1 1 1 1 2 1 1 1 1 3 6 1 1 1 4 7 9 1 1 5 8 10 11 1

上述構造和第3行運算有密切關照

第3行採用矩陣點除操作：

>> tans = 1 2 3 4 5 1 1 6 7 8 1 1 1 9 10 1 1 1 1 11 1 1 1 1 1>> t./tans = 1.0000 0.5000 0.3333 0.2500 0.2000 2.0000 1.0000 0.1667 0.1429 0.1250 3.0000 6.0000 1.0000 0.1111 0.1000 4.0000 7.0000 9.0000 1.0000 0.0909 5.0000 8.0000 10.0000 11.0000 1.0000

分母轉置和第2行索引批量賦值結果矩陣正好前後呼應，銜接天衣無縫。

Pro43011. Convert a vector to a lower triangular matrix

繼續前一個問題的討論，仍然是和矩陣的上、下三角元素有關，變換一個角度：不再要求上三角元素和下三角之間存在關聯，直接賦值為0，但下三角發生一些變化，如題。

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

【原題】：I now have a row vector and I want to convert it to a lower trilangular matrix.

The rows of the lower trilangular matrix have the corresponding elements of the vector.

For example:

input = 1:16output = [ 1 0 0 0 0 0 2 3 0 0 0 0 4 5 6 0 0 0 7 8 9 10 0 0 11 12 13 14 15 0 16 0 0 0 0 0]

【解釋】：前一個問題有3點改動:(1) 主對角線計入輸入數組元素統計；(2). 元素按行排布方式計入；(3). 輸入數組元素可能小於下三角陣元素數量，尾數不足要用0補足。這3處變化足以使得問題的解決方式存在很多潛在的變化。

Test

Code Input

%%x = 1:15;y_correct = [ 1 0 0 0 0 2 3 0 0 0 4 5 6 0 0 7 8 9 10 0 11 12 13 14 15];

assert(isequal(vec2tril(x),y_correct))

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

【分析】：有2個小問題需要解決，首先仍然是根據數組長度判定維度的問題；其次是補0操作何時進行，兩個Test的數組長度正好差1，但結果不同，一個是正好全部填充，一個則由於新增一行，尾數不夠補0。

逐個元素賦值的循環寫法

如下代碼方案使用for循環逐個賦值，這是本問題中，個人感覺寫得比較乾淨整齊的for循環代碼，構造符合條件的全0矩陣採用了floor+sqrt試湊，for循環中，用if保證循環中列坐標小於行坐標時才能用ans(r,c)=i賦值，也就是為下三角元素賦值提供保證。

function ans = vec2tril(x)zeros(floor(sqrt(max(x))) + 2);[r, c] = deal(1);for i = x ans(r, c) = i; c = c + 1; if c > r r = r + 1; c = 1; endend

另一種方式是while循環，判定維度選擇nnz+tril這種比較啰嗦的寫法，原因是我感覺還沒人用過，只是為了保留一個不一樣的命令組合，當然更合理的方法是用cumsum或者其他按照推導公式試湊的表達式：

function ans = vec2tril(x) zeros(find(arrayfun(@(i)nnz(tril(ones(i))),1:10)>=x(end),1)); i=1; while ~isempty(x) t=1:min([numel(x) i]); ans(i,t)=x(t); x(t)=[];i=i+1; endend

批量索引賦值

Alfonso提交了和前一問題類似的，利用索引批量賦值的解法：

function ans = vec2tril(x)[];while length(x) n=min(length(x),length(ans)+1); ans(end+1,1:end+1)=[x(1:n) zeros(1,length(ans)+1-n)]; x(1:n)=;end

while循環第1句「n=min(length(x),length(ans)+1)」是彈性擴維的典型案例，接下來的一句按照索引條件逐行賦值構造；最後一句把輸入x中已賦值元素自x刪除，保證下一次循環能繼續從x(1)開始。

同樣思路，yurenchu給出另一種寫法：

function ans = vec2tril(x) n = find(cumsum(x)>=length(x),1); triu(ones(n)); ans(ans|0) = 1:n*(n+1)/2; ans.*(ans<=max(x))end

選擇cumsum+find的組合計算出輸出矩陣維度，用triu得到上三角矩陣後，轉換為邏輯索引的「ans|0」操作需要mark一下，這句通過邏輯索引賦值的方法之前已經提及，不過有2個問題仍然要提及：首先，邏輯索引賦值順序是列排布，而題目要求的是行排布，顯然需要跟轉置操作；其次，1:n*(n+1)/2是把整個上三角陣元素全部填滿，以Test 2為例，即：x=1:16，第3句賦值得到的結果肯定是超過最大值16的：

>> ansans = 1 2 4 7 11 16 0 3 5 8 12 17 0 0 6 9 13 18 0 0 0 10 14 19 0 0 0 0 15 20 0 0 0 0 0 21

因此後面不僅需要轉置，還需要去掉大於max(x)的元素值，因此最後一句根據轉置矩陣創建0-1同維邏輯矩陣和轉置矩陣本身做點乘。

上述為通過邏輯索引構造賦值，yurenchu還給出了使用低維索引批量賦值的代碼方案：

function ans = vec2tril(x) find(cumsum(x)>=length(x),1); q = find(triu(ones(ans))); zeros(ans); ans(q(x)) = x; ans;end

第2句中的find，內部並沒有邏輯操作，但輸入變數triu(ones(ans))中元素如果大於0，則返回邏輯值1，反之為邏輯值0，這是Cody比較常用的操作。yurenchu這2個代碼全部屬於常見常用操作，不偏不怪，但設計合理、思路明晰，是值得學習的矢量化代碼，所提交的第2段代碼也是該問題目前的leading solution。

劉鵬也用到了索引批量賦值操作，但輸出矩陣維度計算用到的是roots+ceil+max的操作方式。

function y = vec2tril(x) n = ceil(max(roots([1 1 -2*numel(x)]))); y = zeros(n); y(find(triu(true(n)),numel(x))) = x; y = y.;end

MATLAB中提取矩陣元素邏輯索引的操作方式方法多種多樣，例如J.S.Kowontan用到isnan：

function ans = vec2tril(x) triu(NaN(ceil(-0.5+sqrt(0.25+2*length(x))))); ans(find(ans,length(x)))=x; ans(isnan(ans)) = 0; ans;end

稀疏矩陣命令sparse

最後是大神LY Cao的leading solution(和yurenchu以低維索引賦值的解法size相同-並列)：

function ans = vec2tril(x) round(sqrt(2*x)); sparse(ans,cumsum(1-[1 diff(ans)].*(ans-1)),x,ans(end),ans(end));end

短短兩行代碼，但構思如天馬行空毫無斧鑿痕迹，想理解並不容易，先說稀疏矩陣構造命令sparse：

S =sparse(i,j,v,m,n)

選擇5輸入調用格式的sparse命令參數如下：i、j分別代表所構造稀疏矩陣中需要賦非0數值的行列索引號，兩個變數本身必須同維，且和第3參數v維度相同，因為索引位置賦值來自第3參數v，即輸出變數S(i(k),j(k))=v(k)；第4和第5參數m和n則是總輸出變數S的維度，即S是m×n稀疏矩陣。

以x=1:16的輸入，下斷點調出程序運行中間變數，並相應對號入座：

K>> i=ansi = 1 2 2 3 3 3 4 4 4 4 5 5 5 5 5 6K>> j=cumsum(1-[1 diff(ans)].*(ans-1))j = 1 1 2 1 2 3 1 2 3 4 1 2 3 4 5 1K>> v=xv = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16K>> m=ans(end)m = 6

LY Cao和Alfonso兩人令同儕瞠乎而其後的模型構造能力總是令人印象深刻，上述代碼關鍵是i和j兩處的矩陣行列索引構造，他是如何想到用「round(sqrt(2*x))」和「cumsum(1-[1 diff(ans)].*(ans-1))」來構造並取得下三角矩陣的行列索引的？初步猜測大神是由test結果數據反向推斷得到的索引構造代碼，但數據和代碼間的聯繫，我坐這兒把答案平攤著看了半天，也才半懂不懂，自己想出這種解答估計沒指望了。

才知道限制我想像力的除了貧窮原來還有平庸...淚目

延伸討論

下三角矩陣維度計算的總結

最後2道題目中，有一個類型相似的問題，即：通過給定輸入數組長度，即下三角元素個數，判斷輸出矩陣(n×n)的維度，從各路好漢提供的代碼方案發現：解決方案五花八門千奇百怪。當時為湊熱鬧，我也弄了兩個，一個是tril+nnz的組合，用arrayfun獲得1:10維的方陣維度；另一個是根據test輸入的特徵偷雞鑽了個空子。

事實證明，前一個太啰嗦，後一個能偷雞是因為Test只有2個，一個是x=1:15，一個是x=1:16，顯然max(x)-10就是方陣維度，但一方面對其他長度數組不具備適用性，更何況，正確答案的size也未必就能大多少，這種偷雞港真是不咋滴，所以這類方案在下面的總結中就被我選擇性忽略了。

數列求和公式轉化為方程求根

第3個問題的規律很容易看出來，顯然第n項為：

1+2+3+...+n=n(n+1)/2=length(x)

一個為人所熟知的數列求和公式，因此求n就很簡單了，通過求根命令roots對多項式方程：

f(n)=n^2+n-2*length(x)

求根，最大的根向上取整得到的就是所需結果，向上取整是因為輸入變數x的長度小於等於輸出矩陣下三角元素的個數。

n = ceil(max(roots([1 1 -2*numel(x)])))

find命令索引查找

n階方陣下三角元素個數恰好是1+2+3+...+n，這在MATLAB中可以表示為cumsum(1:n)，那麼大於等於length(x)的第1個值的索引，就是希望得到的矩陣維度n，因此：

find(cumsum(x)>=length(x),1)

試湊

無論LY Cao的「round(sqrt(2*x))」，或J.S.Kowontan的「ceil(-0.5+sqrt(0.25+2*length(x)))」運算都體現了數字內在規律在「松」和「緊」之間度的把握。最後其實都形成了各自配合下方構造下三角陣元素構造的特徵數組，結果是一樣的：

>> arrayfun(@(i)ceil(-0.5+sqrt(0.25+2*i)),1:16)ans = 1 2 2 3 3 3 4 4 4 4 5 5 5 5 5 6>> round(sqrt(2*x))ans = 1 2 2 3 3 3 4 4 4 4 5 5 5 5 5 6

正好契合了每行非零元素的個數，而且這個數字的規律，正好是用repelem構造的逆過程：

>> cell2mat(arrayfun(@(i)repelem(i,1,i),1:6,un,0))ans = 1 2 2 3 3 3 4 4 4 4 5 5 5 5 5 6 6 6 6 6 6>> ans(x)ans = 1 2 2 3 3 3 4 4 4 4 5 5 5 5 5 6

只不過實現並不知道構造結果ans到底在第幾項和x的長度吻合，試探起來相對麻煩。

擴維、彈性數組長度以及逗號表達式

關於擴維

第1個題目用到repelem對數組或者元素擴維，MATLAB的數組或矩陣操作中經常會遇到擴維的問題，repelem是MATLAB在2015a中推出的新函數，所以順便聊幾句擴維，擴維的方法非常多，在我和彬彬的書中專門有一節講述了大概7-8種不同的擴維方法，有興趣可以去看看。這裡只是簡單說說repelem、kron這2個函數

Repelem的擴維

這個函數灰強靈活，對第1個參數，按第2個參數給定數量擴展。

>> repelem(1:4,1:4)ans = 1 2 2 3 3 3 4 4 4 4>> repelem(abc,1:3)ans = abbccc

Kron的擴維也是一個道理：

>> kron(ones(1,4),1:4)ans = 1 至 12 列 1 2 3 4 1 2 3 4 1 2 3 4 13 至 16 列 1 2 3 4>> kron(ones(1,3),1:3)ans = 1 2 3 1 2 3 1 2 3>> kron(ones(1,3),randi(10,2))ans = 7 9 7 9 7 9 1 10 1 10 1 10>> kron(ones(2),randi(10,2))ans = 7 8 7 8 8 4 8 4 7 8 7 8 8 4 8 4>> kron(ones(2),abc)ans = 97 98 99 97 98 9997 98 99 97 98 99

這裡需要注意2點：首先，kron是整體擴展，不深入所擴展數組的內部元素，其次，字元串擴展時會自動轉成ascii碼。

關於數組彈性長度

MATLAB中的數組或者矩陣的彈性長度(維度)構造屬於比較靈活的操作，例如避免明顯的判斷流程，用簡潔代碼取得合適的彈性長度向量，這體現了MATLAB易用性，同時也考驗寫代碼時對函數命令理解是否深刻。正如前面幾個問題中對min、max函數的靈活運用就是典型例子。

關於逗號表達式

逗號表達式一個比較常見的用法是結合varargin，做輸入變數數量未定情況下的多變數統一賦值，而不至於出現類似下面的函數開頭：

function [a1, b1,c1,d1,a2,b2,c2,d2,…,an,bn,cn,dn]=Fun(A1, B1,C1,…,An,Bn,Cn,Dn)

這樣的程序簡直叫人人頭昏腦漲，看個開頭，情緒已經有抑鬱的徵兆，改成如下套路顯然好很多：

function v=Test(varargin)[v{1:nargin}]=varargin{:};end

此時輸入的v就和輸入部分有了一定的關聯，例如，調用時可能有如下這種彈性結果：

>> v=Test(1,2)v = 1×2 cell 數組 {[1]} {[2]}>> v=Test(2,randi(10,3),abcd)v = 1×3 cell 數組 {[2]} {3×3 double} {abcd}>> v=Test(2,randi(10,3),abcd,@(x)x.^2)v = 1×4 cell 數組 {[2]} {3×3 double} {abcd} {@(x)x.^2}

這種賦值形式，不是為了酷，因為再酷也不能幫你把妹子，這麼做目的就是大大簡化代碼size，讓輸出變數更易管理和適應不同的輸入狀況。

小結

這篇筆記的內容重點是通過3個問題的多種代碼方案，介紹MATLAB中，數組彈性維度的實現方式，以及根據構造的索引(包括邏輯索引和高、低維索引)批量賦值的具體方法。

在具有彈性長度的數組構造時，應注意使用基本函數例如文中多次提到的min/max等基本命令，一方面可以確保在循環或其他運算中，長度動態改變的數組不至於產生返回空矩陣的錯誤；另一方面這種操作對代碼簡化、避免if判斷起到很大作用。彈性的長度包含2層含義，一方面是數組本身具有變化長度，當然這在MATLAB語言中並不推薦，如果不事先聲明維度而直接對變數賦值，會出現橙色的代碼警告標識，但代碼本身的執行沒有問題；另一方面則是索引具有動態的長度，因此在索引構造和批量賦值中，應考慮構造富有彈性的數組，而在前述3個問題中，多次出現邏輯索引和高、低維索引構造之後的批量賦值，就是這類技巧的具現。