顯示具有 數學 標籤的文章。 顯示所有文章
顯示具有 數學 標籤的文章。 顯示所有文章

2024年10月1日 星期二

完全平方日 20241001

今天山陀兒颱風來襲高雄放假一天, 下午在 MSN 新聞看到這篇 :


原來 20241001 是 4499 的平方啊! 



本世紀上一次是 20151121 (4489 平方), 下一次是 48 年後的 20720704 (4552 平方), 最後一次是 20930625 (4575 平方), 後面這兩個我都遇不到了, 所以山陀兒颱風假的今天是我此生唯一一次注意到的完全平方日啦! 

2023年5月17日 星期三

Python 學習筆記 : 字集排列的過濾 (十七)

本篇為追加過濾規則的最後一個關於字數的限制, 本系列之前的筆記參考 : 


第十五道過濾的規則是 : "每12個字的排列中,這12字加起來的字母總數介於64~85個字母". 此規則不需要用到正規式, 用 Python 的 len() 函式就能搞定. 

第一個測試語料是字母總數少於 64 個字元的情況 (應排除) :

#  perm_test_15.py
#  字母總數須介於 64~85 個字母
words=('air', 'root', 'arch', 'try', 'close', 'amber', 'silent', 'beauty', 'trial', 'engage', 'more', 'after')
chars=''.join(words)
print(chars)
print(len(chars))
if len(chars) < 64 or len(chars) > 85:  # 須 64~85 個字母才入選                      
    print('the permutation is excluded')       
else:      
    print('the permutation is included')     

執行結果為排除 (正確) :

>>> %Run perm_test_15.py   
airrootarchtrycloseambersilentbeautytrialengagemoreafter
56
the permutation is excluded

第二個測試語料是字母總數等於 64 個字元的情況 (應入選) :

#  perm_test_15.py
#  字母總數須介於 64~85 個字母
words=('air', 'root', 'machine', 'jacket', 'close', 'country', 'silent', 'beauty', 'trial', 'engage', 'more', 'after')
chars=''.join(words)
print(chars)
print(len(chars))
if len(chars) < 64 or len(chars) > 85:  # 須 64~85 個字母才入選                      
    print('the permutation is excluded')       
else:      
    print('the permutation is included')    

執行結果為入選 (正確) :

>>> %Run perm_test_15.py   
airrootmachinejacketclosecountrysilentbeautytrialengagemoreafter
64
the permutation is included

第三個測試語料是字母總數等於 85 個字元的情況 (應入選) :

#  perm_test_15.py
#  字母總數須介於 64~85 個字母
words=('abandon', 'negative', 'machine', 'jacket', 'champion', 'country', 'balance', 'beauty', 'obscure', 'engage', 'shoulder', 'identify')
chars=''.join(words)
print(chars)
print(len(chars))
if len(chars) < 64 or len(chars) > 85:  # 須 64~85 個字母才入選                      
    print('the permutation is excluded')       
else:      
    print('the permutation is included')
   
執行結果為入選 (正確) :

>>> %Run perm_test_15.py   
abandonnegativemachinejacketchampioncountrybalancebeautyobscureengageshoulderidentify
85
the permutation is included

第四個測試語料是字母總數多於 85 個字元的情況 (應排除) :

#  perm_test_15.py
#  字母總數須介於 64~85 個字母
words=('abandon', 'negative', 'machine', 'jacket', 'champion', 'country', 'balance', 'kangaroo', 'obscure', 'engage', 'shoulder', 'identify')
chars=''.join(words)
print(chars)
print(len(chars))
if len(chars) < 64 or len(chars) > 85:  # 須 64~85 個字母才入選                      
    print('the permutation is excluded')       
else:      
    print('the permutation is included')
   
執行結果為排除 (正確) :

>>> %Run perm_test_15.py   
abandonnegativemachinejacketchampioncountrybalancekangarooobscureengageshoulderidentify
87
the permutation is excluded


2023-05-19 補充 :

規則改為這 12 個字用空格串起來後長度要在 64~88 個之間, 故上面的測試程式改成如下 : 

chars=' '.join(words)     # 以空格串接排列中的 12 個字母
print(chars)
print(len(chars))
if len(chars) < 64 or len(chars) > 88:  # 含空格長度須 64~88 個字母才入選                      
    print('the permutation is excluded')       
else:      
    print('the permutation is included')

2023年5月13日 星期六

Python 學習筆記 : 字集排列的過濾 (十)

本篇為追加過濾兩個字尾母音規則的第二個測試, 本系列之前的筆記參考 : 


第八道過濾的規則是 : "每12字排列中,a, i, o, u 四個母音字母(變成一群),以這四個母音字母「結尾」的字有0~1個". 此規則的正規式如下 :

ptn=re.compile(r'\b\w*[aiou]\b')      # 以 a, i, o, u 四個母音字母結尾的字     

第一個測試語料是沒有以 a, i, o, u 四個母音字母結尾的字的情況 (應入選) :

>>> import re 
>>> words=('please', 'age', 'pole', 'check', 'close', 'open', 'silent', 'apology', 'trial', 'engage', 'more', 'license')   
>>> ptn=re.compile(r'\b\w*[aiou]\b')      # a, i, o, u 四個母音字母結尾的字
>>> end_aiou=[True if re.match(ptn, w) else False for w in words]   
>>> end_aiou    
[False, False, False, False, False, False, False, False, False, False, False, False]
>>> if end_aiou.count(True) > 1:    # True 個數須不超過 1 個 才入選                      
    print('the permutation is excluded')       
else:      
    print('the permutation is included') 
    
the permutation is included  

可見含有 0 個以 a, i, o, u 四個母音字母結尾的字被納入. 

其次測試含有 1 個以 a, i, o, u 四個母音字母結尾的語料 : 

>>> words=('please', 'ago', 'pole', 'check', 'close', 'open', 'silent', 'apology', 'trial', 'engage', 'more', 'license')   
>>> ptn=re.compile(r'\b\w*[aiou]\b')      # a, i, o, u 四個母音字母結尾的字
>>> end_aiou=[True if re.match(ptn, w) else False for w in words]   
>>> end_aiou    
[False, True, False, False, False, False, False, False, False, False, False, False]
>>> if end_aiou.count(True) > 1:    # True 個數須不超過 1 個 才入選                      
    print('the permutation is excluded')       
else:      
    print('the permutation is included') 
    
the permutation is included 

可見含有 1 個以 a, i, o, u 四個母音字母結尾的字也被納入. 

接下來測試含有 2 個以 a, i, o, u 四個母音字母結尾的語料 : 

>>> words=('please', 'ago', 'pole', 'check', 'close', 'open', 'silent', 'apology', 'trial', 'engage', 'more', 'alumni')   
>>> ptn=re.compile(r'\b\w*[aiou]\b')      # a, i, o, u 四個母音字母結尾的字
>>> end_aiou=[True if re.match(ptn, w) else False for w in words]   
>>> end_aiou    
[False, True, False, False, False, False, False, False, False, False, False, True]
>>> if end_aiou.count(True) > 1:    # True 個數須不超過 1 個 才入選                      
    print('the permutation is excluded')       
else:      
    print('the permutation is included') 
    
the permutation is excluded    

以上測試可以寫成如下程式於命令列執行 :

#  perm_test_8.py
#  以 a, i, o, u 四個母音字母結尾的字不超過 1 個
import re

words=('please', 'age', 'pole', 'check', 'close', 'open', 'silent', 'apology', 'trial', 'engage', 'more', 'license')
ptn=re.compile(r'\b\w*[aiou]\b')      # a, i, o, u 四個母音字母結尾的字
end_aiou=[True if re.match(ptn, w) else False for w in words]
print(end_aiou)
print(end_aiou.count(True))
if end_aiou.count(True) > 1:    # True 個數須不超過 1 個 才入選                      
    print('the permutation is excluded')       
else:      
    print('the permutation is included')     

2023年5月8日 星期一

Python 學習筆記 : 字集排列的過濾 (八)

這兩天重新審視了字集過濾的規則, 修正了規則 5 以免規則 6 形同具文, 為了驗證這些規則的正確性, 使用有限的 14 個排列來測試, 以下是測試紀錄. 本系列之前的文章參考 :

本系列之前的文章參考 : 

測試用的 14 個排列語料如下 (檔名 permutation_test_in.csv) : 

ability,give,about,echo,ugly,idea,fake,grow,close,brick,merge,nature
ability,give,about,echo,ugly,idea,fake,grow,close,obey,merge,nature
sugar,frown,pole,million,hair,close,silent,apology,engage,dish,harvest,license
abandon,ability,bag,balance,key,scan,select,shock,radar,radio,ugly,umbrella
abandon,baby,cabbage,dad,eager,fabric,gadget,habit,ice,jacket,kangaroo,lab
sugar,frown,pole,million,hair,close,silent,apology,engage,dish,harvest,license
cinnamon,merge,more,memory,grow,anchor,auto,major,push,desk,mass,swallow
knife,brick,quote,interest,kind,jealous,afraid,jar,job,much,hat,umbrella
cost,brick,gate,interest,depart,jealous,afraid,diamond,merit,much,hat,umbrella
cost,brick,gate,interest,depart,face,afraid,diamond,merit,much,hat,umbrella
give,fade,pole,million,gift,close,silent,apology,engage,face,ginger,license
give,fade,pole,million,date,close,silent,apology,engage,face,ginger,license
cheap,fade,pole,check,cheese,close,silent,apology,engage,this,thing,write
cheap,fade,pole,mail,cheese,close,silent,apology,engage,this,thing,write

測試程式如下 (檔名 words_permutation_test.py) :

import re
import time

start=time.time()
# 資料前處理 : 讀取 CSV 檔轉成串列
with open('permutation_test_in.csv', 'r', encoding='utf8') as fr:
    with open('permutation_test_out.csv', 'w', encoding='utf8') as fw:
        lines=fr.readlines()
        i=1 # 排列計數器
        for line in lines:
            words=line.replace('\n', '')
            words=words.split(',')
            print(i, end=':')  # 印出排列數
            print(words)       # 印出排列 (12 字的 tuple)
            i=i+1  # 排列數增量 1
            #rule1 : 母音字母 (a, e, i, o, u) 開頭的字最多出現 5 次
            ptn=re.compile('^[aeiou].*') # 母音字母開始
            first=[w[0] for w in words if re.match(ptn, w)]
            if len(first) > 5: # 母音開頭字超過 5 次
                print(" : rule 1 excluded")
                continue
            #rule2 : 相同字母開頭的字最少 1 組, 最多 4 組, 母音與子音可同時併計
            first=[w[0] for w in words]  # 找出各字之開頭字母串列
            first_diff=list(set(first))  # 找出不同開頭字母串列
            fc=[first.count(fd) for fd in first_diff]
            fc1=[True if c > 1 else False for c in fc]
            if fc1.count(True) < 2 or fc1.count(True) > 4: 
                print(" : rule 2 excluded")
                continue
            # rule3: 相同字母開頭的字最少 2 個, 最多 4 個
            first=[w[0] for w in words]
            first_diff=list(set(first))
            fc=[first.count(fd) for fd in first_diff]
            fc1=[True if first.count(fd) > 1 else False for fd in first_diff]
            fc2=[True if first.count(fd) > 4 else False for fd in first_diff]
            if fc1.count(True) < 1 or fc2.count(True) > 0: 
                print(" : rule 3 excluded")
                continue  
            # rule4: 以 j, k, q, y, z 開頭的字最多只能有 1 個
            ptn=re.compile('^[jkqyz].*')
            first=[w[0] for w in words  if re.match(ptn, w)]
            first_diff=list(set(first))
            fc=[first.count(fd) for fd in first_diff]
            if sum(fc) >= 2:         
                print(" : rule 4 excluded")
                continue
            # rule5: 前 2 個字母開頭相同的字最多只能出現 2 次 (sh, ch, th, wr, un 例外)
            ptn='^(?!(sh|ch|th|wr|un))[a-zA-Z]*'
            first2=[w[0:2] for w in words if re.match(ptn, w)]
            first2_diff=list(set(first2))
            f2c=[True if first2.count(fd) > 2 else False for fd in first2_diff]
            if f2c.count(True) > 0: # 
                print(" : rule 5 excluded")
                continue
            # rule6: 以 sh, ch, th, wr, un 開頭的字, 前三個字母相同者不能超過 2 個
            ptn=re.compile('^(sh|ch|th|wr|un)[a-zA-Z]*')
            first3=[w[0:3] for w in words if re.match(ptn, w)]
            first3_diff=list(set(first3))
            f3c=[True if first3.count(fd) > 2 else False for fd in first3_diff]
            if f3c.count(True) > 0: # 
                print(" : rule 6 excluded")
                continue
            # 通過上面 6 個過濾 : 存入檔案
            str=' '.join(words)
            print(str)
            fw.write(str + '\n')        
end=time.time()
print(f'time elapsed : {end-start}')

執行結果如下 :

>>> %Run words_permutation_test.py   
1:['ability', 'give', 'about', 'echo', 'ugly', 'idea', 'fake', 'grow', 'close', 'brick', 'merge', 'nature']
ability give about echo ugly idea fake grow close brick merge nature
2:['ability', 'give', 'about', 'echo', 'ugly', 'idea', 'fake', 'grow', 'close', 'obey', 'merge', 'nature']
 : rule 1 excluded   
3:['sugar', 'frown', 'pole', 'million', 'hair', 'close', 'silent', 'apology', 'engage', 'dish', 'harvest', 'license']
sugar frown pole million hair close silent apology engage dish harvest license
4:['abandon', 'ability', 'bag', 'balance', 'key', 'scan', 'select', 'shock', 'radar', 'radio', 'ugly', 'umbrella']
 : rule 2 excluded   
5:['abandon', 'baby', 'cabbage', 'dad', 'eager', 'fabric', 'gadget', 'habit', 'ice', 'jacket', 'kangaroo', 'lab']
 : rule 2 excluded   
6:['sugar', 'frown', 'pole', 'million', 'hair', 'close', 'silent', 'apology', 'engage', 'dish', 'harvest', 'license']
sugar frown pole million hair close silent apology engage dish harvest license
7:['cinnamon', 'merge', 'more', 'memory', 'grow', 'anchor', 'auto', 'major', 'push', 'desk', 'mass', 'swallow']
 : rule 3 excluded    
8:['knife', 'brick', 'quote', 'interest', 'kind', 'jealous', 'afraid', 'jar', 'job', 'much', 'hat', 'umbrella']
 : rule 4 excluded   
9:['cost', 'brick', 'gate', 'interest', 'depart', 'jealous', 'afraid', 'diamond', 'merit', 'much', 'hat', 'umbrella']
cost brick gate interest depart jealous afraid diamond merit much hat umbrella
10:['cost', 'brick', 'gate', 'interest', 'depart', 'face', 'afraid', 'diamond', 'merit', 'much', 'hat', 'umbrella']
cost brick gate interest depart face afraid diamond merit much hat umbrella
11:['give', 'fade', 'pole', 'million', 'gift', 'close', 'silent', 'apology', 'engage', 'face', 'ginger', 'license']
 : rule 5 excluded   
12:['give', 'fade', 'pole', 'million', 'date', 'close', 'silent', 'apology', 'engage', 'face', 'ginger', 'license']
give fade pole million date close silent apology engage face ginger license
13:['cheap', 'fade', 'pole', 'check', 'cheese', 'close', 'silent', 'apology', 'engage', 'this', 'thing', 'write']
 : rule 6 excluded
14:['cheap', 'fade', 'pole', 'mail', 'cheese', 'close', 'silent', 'apology', 'engage', 'this', 'thing', 'write']
cheap fade pole mail cheese close silent apology engage this thing write
time elapsed : 0.11959671974182129

可見規則 1~6 都有正確過濾 (7 個被濾掉), 輸出檔內容如下 :

ability give about echo ugly idea fake grow close brick merge nature
sugar frown pole million hair close silent apology engage dish harvest license
sugar frown pole million hair close silent apology engage dish harvest license
cost brick gate interest depart jealous afraid diamond merit much hat umbrella
cost brick gate interest depart face afraid diamond merit much hat umbrella
give fade pole million date close silent apology engage face ginger license
cheap fade pole mail cheese close silent apology engage this thing write

14 個排列有 7 個被規則 1~6 過濾掉, 剩下 7 個. 

2023年4月17日 星期一

Python 學習筆記 : 用 itertools.permutations() 排列字集

以前高師大英語所碩班的同學明中昨晚來電詢問 : 如果有 1500 個字的英文字集每次挑 12 個字出來排列會有多少種不同的排列? 電腦的算力能否處理? 關於第一個問題, Python 的內建模組 math 內有一個 perm() 函式可輕易計算出來 (> v3.8) : 

Python 3.11.2 (C:\Users\User\AppData\Local\Programs\Python\Python311\python.exe)
>>> import math   
>>> math.perm(1500, 12)     
124147257394529035596269620244764800000
>>> print("%e" %math.perm(1500, 12))       # 用科學表示法呈現
1.241473e+38

哇, 是 38 次方等級的龐然大物, 不知 PC 能否扛得起哩. 不過還會加上一些限制條件, 或許不會到這麼龐大啦 (其實還是很大), 參考 :


此處摘要一下數學中排列 (permutation) 與組合 (combination) 的差別, 排列與順序有關, 組合則與順序無關, 元素一樣但順序不一樣算不同的排列, 但卻是同一種組合, 例如 a, b, c 三個字元每次挑兩個, 有 ab, ba, ac, ca, bc, cb 六種排列, 但卻只有 ab, ac, bc 三種組合 (因為 ab/ba, ac/ca, bc/cb 都各算一種組合), 其數學公式如下 : 



雖然可以用 math.factorial() 函式根據上面的算式來計算排列組合數, 但在 Python 3.8 版後 math 模組新增了 perm()comb() 函式, 可用來分別計算排列數與組合數, 不需要自行套公式來算, 參考 :


組合數因為與順序無關會比排列數來得少, 1500 取 12 結果是 29 次方等級 :

>>> math.comb(1500, 12)   
259179212333589356687471649875   
>>> print("%e" %math.comb(1500, 12))    
2.591792e+29

在 1500 個字集準備好之前, 我先用 A~Z 這 26 個字母做個簡單測試, 作法是利用 Python 標準函式庫 itertools 中的 permutations() 函式來列出所有排列, 參考 : 


itertools.permutations() 函式可傳入兩個參數 : 


第一個參數為可迭代物件 (必要參數), 第二個參數 (備選參數) 是傳回值的排列長度 (即選幾個來排列), 例如從 26 個英文字母中取 5 個來排列 : 

>>> import math  
>>> from itertools import permutations     
>>> letters=list('ABCDEFGHIJKLMNOPQRSTUVWXYZ')     
>>> letters   
['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z']
>>> math.perm(26, 5)      # 26 取 5 有近 79 萬種排列方式
7893600
>>> for p in permutations(letters, 5): 
    print(p)                          # 會印很久
    
將近 79 萬種排列在互動環境就跑了好久啊! 真不敢想像那 38 次方會怎樣! 不過, 在互動環境 print 結果本來就慢很多, 我將上面的指令加上計時功能寫成如下程式檔案, 並將排列結果輸出到檔案, 在 Thonny 中執行就快多了, 也才用了 30 秒不到 : 

# perm_test.py
import time
import math
from itertools import permutations

start=time.time()
letters=list('ABCDEFGHIJKLMNOPQRSTUVWXYZ')
f=open('permutation.csv', 'w')
for p in permutations(letters, 5):
    f.write(','.join(p) + '\n')
f.close()
end=time.time()
print(f'time elapsed : {end-start}')

>>> %Run perm_test.py
time elapsed : 27.593969345092773   

結果得到一個大小約 84MB 的 csv 檔, 節錄內容頭尾如下 :


.......


雖然看來很快, 但還無法估算 38 次方要跑多久才跑得完, 以及檔案有多大. 


2023-04-20 補充 :

昨天收到測試用的 238 字字集, 存在 A~Z 欄的 Excel 檔案裡, 我將其另存成 CSV 檔 : 




其中連續的逗號來自空格, 清除換列字元與空格後, 實際字數為 238 字. 如果從 238 個字中挑選 12 個, 其排列數計算如下 :

>>> import math   
>>> math.perm(238, 12)    
24917297304498261278500684800
>>> print("%e" %math.perm(238, 12))    
2.491730e+28

哇, 是 18 次方等級的. 

接著寫了如下之程式來處理 (238, 12) 的排列情況 :

# words_permutation.py
import time
import math
from itertools import permutations

start=time.time()
with open('TW_seed_phrase_test.csv', 'r', encoding='utf8') as f:
    lines=f.readlines()
words_str=','.join(lines)
words_str=words_str.replace('\n', '')
words_list=words_str.split(',')
words_list=list(filter(lambda word: word != '', words_list))
words_list.sort()
with open('words_permutation.csv', 'w', encoding='utf8') as f:
    for p in permutations(words_list, 12):
        f.write(' '.join(p) + '\n')
print(words_list)
print(len(words_list))
end=time.time()
print(f'time elapsed : {end-start}')

此處先用字串的 split() 方法以逗號為界將字串 words_str 拆成串列, 但 EXCEL 中的空格會在匯出為 CSV 檔時產生連續逗號, 拆解時會在串列中變成空字串, 這裡利用了內建函式 filter() 來清除這些空字串, 它的傳入參數是一個傳回 True/Fasle 的函式 (辨認可迭代物件之元素是否為空字串), filter() 會將傳回 True 的可迭代物件元素傳回, 參考 : 


然後將請裡過的 words_list 串列丟給 math.permutations() 去做排列, 它會傳回一個 generator, 只要迭代此 generator 就可以產生各個排列結果 (tuple), 將此 tuple 以空格隔開組成字串後寫入輸出檔案即可. 關於檔案讀寫參考 :


跑了 5 個多小時還沒結束 ...... 持續觀察中. 

跑了快 6 小時發現程式因為 "maximum recursion depth exceeded" (超過最大迭代深度) 跳離, 同時 400 多 GB 的硬碟也快用光了 : 

Python 3.11.2 (C:\Users\User\AppData\Local\Programs\Python\Python311\python.exe)
>>> object address  : 00000229B953BDC0
object refcount : 2
object type     : 00007FF838D13550
object type name: RecursionError 
object repr     : RecursionError('maximum recursion depth exceeded')
lost sys.stderr

Process ended with exit code 1.

檢視輸出檔真的來到 424GB 左右 : 




用 EXCEL 開啟時會詢問 CSV 的區隔字元, 可看出輸出結果正確, 只是檔案太大開了好久都沒開起來只好放棄 : 




看樣子必須想辦法切割輸出的 CSV 檔了.

2021年12月3日 星期五

好書 : 數學, 確定性的失落

今天在研究 Python 的 type 時找到下面這篇文章 : 


作者是學數學的, 他談到了數學的公理系統 (公理 axioms 是雖然不知如何證明為真, 但經由人類長久實踐的經驗, 被普遍認為是不證自明的基本命題), 推薦我們去讀一讀美國數學史學家, 紐約大學教授 Morris Kline 寫的 "數學, 確定性的失落 (Mathematics: The Loss of Certainty)" 這本書 : 


Source : 博客來


雖然不確定自己是否能看得下去 (確定性的失落 XD), 但這應該是一本好書. 市圖有收藏這本, 所以不用買, 那就借來看看唄. 

2021年10月25日 星期一

二哥的機率與統計課本

二哥說這學期機率與統計有點難, 我查了他傳給我的書單, 老師用的教科書是多倫多大學教授 Alberto Leon-Garcia 於 2009 年出版的 "機率統計與隨機過程" 第三版 (Prentice Hall) : 



Source : 天瓏


此書是專為電機工程應用而寫, 故書中範例以電機與資訊工程相關問題為主, 此版還添加了 MATLAB/Octave 範例程式, 雖然沒有很多. 第三版最主要的特色是將離散隨機變數與連續隨機變數拆開在三四章, 同時新增第八章介紹統計學, 第十一章則完整探討離散與連續馬可夫鏈. 

感覺能出到第三版的教科書應該都還不錯, 我以前工程統計學得不怎樣, 拿這本來複習好像不錯, 但若書中範例能用 Python 更好 (有想要邊讀邊給它改寫範例程式, 這又得去了解 Matlab 與 Python 的對應才行).  

2021年10月21日 星期四

Python 學習筆記 : SciPy 的線性代數子套件 linalg 測試

SciPy 的線性代數子套件 linalg 是以 Numpy 的 ndarry 陣列為基礎建構的, 它提供了比 Numpy 更多更完整的線性代數函式庫, 事實上 SciPy 可視為 Numpy 的擴充集 (superset), 但 Numpy 基於歷史緣故仍保有核心的線性代數函式. 


參考書籍 :

基礎線性代數 (黃學亮, 五南, 2013)

使用 SciPy 線性代數子套件 linalg 的方式之一是匯入所有 linalg 的函數 : 

from scipy.linalg import *   

此方式很方便, 因為可以直接呼叫此子套件的函數, 但要注意命名空間汙染問題. 

比較安全的方式是 :

from scipy import linalg  

然後用 linalg 去呼叫旗下的函數, 例如 linalg.inv(A).  

常用函數如下表 : 


 scipy.linalg 常用函數 說明
 det(A) 傳回方陣 A 的行列式值
 inv(A) 傳回方陣 A 的反矩陣
 solve(A, b) 傳回線性方程組 Ax=b 的變數向量 x 之解 (串列)
 eig(A) 傳回方陣 A 的特徵值 (eigen value) 與特徵向量 (eigen vector)
 qr(A) 將方陣 A 進行 QR 分解, 傳回正交矩陣 Q 與上三角矩陣 (元組)
 lu(A) 將方陣分解為置換矩陣 P, 上三角矩陣 U, 下三角矩陣 L (元組)


參考 scipy.linalg 子套件的說明文件 : 



1. 用 det() 函數求方陣的行列式 (determinant) 值 : 
 
一個二階方陣 A 的行列式的算法為 :


即左上右下對角線元素相乘後相加, 再減所有右上左下對角線元素之乘積. 同樣的算法在三階方陣 A 的行列式為 :




下面用 linalg 的 det() 函數來計算上面兩個方陣的行列式 :

>>> A=np.array([[1, -1],[1, 2]])      # 二階方陣
>>> linalg.det(A)        # 計算方陣 A 的行列式值
3.0
>>> A=np.array([[3, 2, 3], [2, 5, -7], [1, 2, -2]])     # 三階方陣
>>> linalg.det(A)        # 計算方陣 A 的行列式值
3.0000000000000013   

下面這個三階方陣的行列式值為 0, 例如 : 
 
>>> A=np.array([[2, 2, 2], [2, 2, 2], [2, 2, 2]])     # 三階方陣
>>> A   
array([[2, 2, 2],
       [2, 2, 2],
       [2, 2, 2]])
>>> linalg.det(A)    
0.0

矩陣的行列式有如下特性 :
  • 行列式值為 0 的矩陣為不可逆, 亦即其反矩陣不存在. 
  • 若方陣 A 任何一列或一行均為 0, 則 det(A)=0.
  • 若方陣 A 任何兩列值相同, 則 det(A)=0.
  • 若方陣 A 為對角, 上三角或下三角矩陣, 則 det(A)=對角線元素的乘積.
  • 任何階的單位矩陣其行列式值必為 0.
  • 兩個方陣內積的行列式值等於個別方陣行列式之乘積 det(AB)=det(A)*det(B).
  • 方陣 A 轉置後的行列式值等於方陣 A 的行列式值. 
  • 方陣 A 的 n 次方 (n 次內積) 的行列式值為 A 行列式值的 n 次方. 


2. 用 inv() 求方陣的反矩陣 (inverse matrix) : 

在線性代數中, 一個 n 階方陣 A 的反矩陣 B 的定義是, A 與 B 的內積為 n 階單元方陣 :  




A 的反矩陣 B 記做 :


一個方陣的反矩陣若存在, 此矩陣稱為可逆 (invertible), 可逆的矩陣又稱為非奇異矩陣 (non-singular matrix); 反之, 若方陣的反矩陣不存在即為不可逆, 不可逆的方陣被稱為奇異矩陣 (singular matrix). 反矩陣若存在, 則它是唯一的, 即一個方陣不可能有兩個反矩陣. 

線性代數中反矩陣是透過伴隨矩陣 (adjugate matrix) 計算出來的, 反矩陣與伴隨矩陣存在一個比例關係, 將方陣 A 的伴隨矩陣 A* 除以 A 的行列式 |A| 即可得到反矩陣, 因此行列式值為 0 之矩陣無法求反矩陣 : 


SciPy 子套件 linalg 的 inv() 函數可用來求方陣的反矩陣, 例如 : 

>>> import numpy as np    
>>> from scipy import linalg    
>>> A=np.array([[1, -1],[1, 2]])          # 定義一個方陣 A
>>> A   
array([[ 1, -1],
       [ 1,  2]])
>>> B=linalg.inv(A)           # 求 A 的反矩陣 B
>>> B   
array([[ 0.66666667,  0.33333333],
       [-0.33333333,  0.33333333]])
>>> np.dot(A, B)                                              # 計算 A, B 的內積
array([[1.00000000e+00, 0.00000000e+00],
       [1.11022302e-16, 1.00000000e+00]])
>>> np.dot(B, A)                                              # 計算 B, A 的內積
array([[ 1.00000000e+00, -1.11022302e-16],
       [ 0.00000000e+00,  1.00000000e+00]])    

A, B 交換做內積結果都是單元矩陣, 可見 A 的反矩陣為 B. 

反矩陣的用途之一是用來解線性聯立方程組 AX=b, 其中 A 為係數矩陣 (coefficient matrix), X 是變數向量, b 則是右手係數 (right-hand coefficients) 向量, 則方程組的解如下 : 


例如下面這個三元線性聯立方程組 : 


寫成矩陣方程式為 :


此方程組的解為係數矩陣之反矩陣與右手係數向量之內積 :


下面用 linalg 的 inv() 函數來計算其解 : 

>>> A=np.array([[3, 2, 3], [2, 5, -7], [1, 2, -2]])       # 係數矩陣
>>> A   
array([[ 3,  2,  3],
       [ 2,  5, -7],
       [ 1,  2, -2]])
>>> b=np.array([9, -12, -3])         # 右手係數向量
>>> b   
array([  9, -12,  -3])
>>> B=linalg.inv(A)       # 係數矩陣的反矩陣
>>> B 
array([[ 1.33333333,  3.33333333, -9.66666667],
       [-1.        , -3.        ,  9.        ],
       [-0.33333333, -1.33333333,  3.66666667]])
>>> np.dot(B, b)             # 聯立方程組的解
array([ 1.00000000e+00, -7.10542736e-15,  2.00000000e+00])

第二個變數之解近乎 0, 故其解為 [x, y, z]=[1, 0, 2]. 

注意, 行列式值為 0 的矩陣不可逆, 沒有反矩陣, 例如 :

>>> A=np.array([[2, 2, 2], [2, 2, 2], [2, 2, 2]])    
>>> A  
array([[2, 2, 2],
       [2, 2, 2],
       [2, 2, 2]])
>>> linalg.det(A)        # 行列式值為 0 的矩陣不可逆
0.0   
>>> linalg.inv(A)       
Traceback (most recent call last):
  File "<pyshell>", line 1, in <module>
  File "C:\Python37\lib\site-packages\scipy\linalg\basic.py", line 979, in inv
    raise LinAlgError("singular matrix")
numpy.linalg.LinAlgError: singular matrix   

行列式值為 0 的矩陣為奇異矩陣, 沒有反矩陣. 


3. 用 solve() 函數求解線性聯立方程組 : 

上面我們使用反矩陣求解線性聯立方程組之解, 事實上呼叫 solve(A, b) 函數並將係數矩陣 A 與右手係數向量 b 傳入分別做為第一與第二參數亦可求解, 它會傳回以串列表示的變數向量, 例如 :

>>> A=np.array([[3, 2, 3], [2, 5, -7], [1, 2, -2]])       # 係數矩陣
>>> A   
array([[ 3,  2,  3],
       [ 2,  5, -7],
       [ 1,  2, -2]])
>>> b=np.array([9, -12, -3])         # 右手係數向量
>>> b   
array([  9, -12,  -3])   
>>> x=linalg.solve(A, b)      # 求解線性聯立方程組
>>> x    
array([1., 0., 2.])        # 變數向量 [x, y, z] 之解

可見答案與上面用反矩陣計算的相同. 


4. 用 eig() 函數求特徵值與特徵向量 :

一個 n 階方陣 A 的特徵值與特徵向量定義如下 : 存在一個非 0 向量 X 與純量 λ 使得 AX=λX, 則 λ 稱為方陣 A 的特徵值, 而向量 X 為其特徵向量. AX=λX 可看成向量 X 經過 A 的線性轉換後得到一個新向量 Y=AX, 而此新向量與原向量 X 呈平行關係, 特徵值即兩平行向量之比例係數. 

AX=λX 經移項後為 (A-λI)X=0, 故 λ 為 A 的特徵值之充要條件為 |A-λI|=0, 將此式展開即可得到 λ 的特徵方程式, 解此方程式之根即可得到特徵值 λ. 呼叫 linalg 子套件的 eig(A) 函數會傳回兩個元素的 tuple, 第一個傳回值為特徵值 λ, 第二個傳回值為特徵向量 : 

l, v=linalg.eig(A)  

例如 :

>>> A=np.array([[1, 2],[3, 2]])     # 二階方陣
>>> A   
array([[1, 2],
       [3, 2]])
>>> l, v=linalg.eig(A)     # 求特徵值與特徵向量
>>> l                                # 有兩個特徵值 : -1 與 4
array([-1.+0.j,  4.+0.j])
>>> v                               # 有兩個特徵向量 : 前者對應特徵值 -1, 後者對應 4
array([[-0.70710678, -0.5547002 ],
       [ 0.70710678, -0.83205029]])

此例有兩組特徵向量, 分別對應兩個特徵值.

特徵方程式可能會有複數根, 例如 : 

>>> A=np.array([[3, 2, 3], [2, 5, -7], [1, 2, -2]])    
>>> A   
array([[ 3,  2,  3],   
       [ 2,  5, -7],
       [ 1,  2, -2]])
>>> l, v=linalg.eig(A)    
>>> l                            # 有三個特徵值
array([4.90057187+0.j        , 0.54971406+0.55676557j,
       0.54971406-0.55676557j])
>>> v                           # 有三組特徵向量 
array([[-0.86169137+0.j        , -0.74295072+0.j        ,
        -0.74295072-0.j        ],
       [-0.44017844+0.j        ,  0.63299867-0.06720996j,
         0.63299867+0.06720996j],
       [-0.25244984+0.j        ,  0.18481478-0.09307649j,
         0.18481478+0.09307649j]])

此例有三個特徵向量, 分別對應三個特徵值. 


5. 用 qr() 函數對方陣進行 QR 分解 :

一個實數方陣 A 的 QR 分解是指它可拆解為一個正交矩陣 Q 與一個上三角矩陣 R 的內積, 所謂正交矩陣意指轉置矩陣等於反矩陣, 亦即正交矩陣與其轉置矩陣之內積為單元矩陣 : 


參考 :


linalg 子套件的 qr(A) 函數會傳回方陣 A 的正交矩陣 Q 與一個上三角矩陣 R 組成之 tuple, 使得 A=QR, 例如 : 

>>> A=np.array([[1, 2],[3, 2]])      # 二階方陣
>>> A   
array([[1, 2],
       [3, 2]])
>>> Q, R=linalg.qr(A)         # 對方陣 A 進行 QR 分解
>>> Q                                    # 正交矩陣 Q
array([[-0.31622777, -0.9486833 ],
       [-0.9486833 ,  0.31622777]])
>>> R                                    # 上三角矩陣 R
array([[-3.16227766, -2.52982213],
       [ 0.        , -1.26491106]])
>>> np.dot(Q, Q.T)              # 正交矩陣與其轉置矩陣之內積必為單元矩陣 I
array([[ 1.0000000e+00, -5.8339229e-18],
       [-5.8339229e-18,  1.0000000e+00]])
>>> B=linalg.inv(Q)             # 求正交矩陣 Q 的反矩陣 B
>>> B                                     # 正交矩陣的反矩陣即為本身之轉置
array([[-0.31622777, -0.9486833 ],
       [-0.9486833 ,  0.31622777]])
>>> np.allclose(A, np.dot(Q, R))      # 檢查 A 與 QR 是否相等
True  

由上面驗算可知, 經 QR 分解後得到的正交矩陣與其轉置矩陣之內積為單元矩陣 I. 此外也使用了 Numpy 的 allclose() 函數檢驗 A 確實等於 QR. 函數 np.allclose() 常用來檢驗兩個陣列 A, B 的每個對應元素是否都相等 (element-wise equal), 若全部相等傳回 True, 否則傳回 False :

np.allclose(A, B) 

其預設判別條件如下 : 


其中 a, b 分別是矩陣 A, B 中的相對應元素, atol 為絕對容忍參數, rtol 為相對容忍參數 (相對於 b 的絕對值). atol 預設值為 1e-8 (即 0.00000001), rtol 預設值為 1e-5 (即 0.00001). 當對應元素 a, b 差之絕對值小於等於右方條件值時傳回 True, 否則傳回 False, 參考 : 


例如 :

>>> np.allclose([1e10,1e-7], [1.00001e10,1e-8])    
False
>>> np.allclose([1e10,1e-8], [1.00001e10,1e-9])   
True

下面是三階方陣的範例 : 

>>> A=np.array([[3, 2, 3], [2, 5, -7], [1, 2, -2]])   # 三階方陣
>>> A   
array([[ 3,  2,  3],
       [ 2,  5, -7],
       [ 1,  2, -2]])
>>> Q, R=linalg.qr(A)        # 對方陣 A 進行 QR 分解
>>> Q                                    # 正交矩陣 Q
array([[-0.80178373,  0.59152048, -0.08512565],
       [-0.53452248, -0.77352678, -0.34050261],
       [-0.26726124, -0.22750788,  0.93638218]])
>>> R                                    # 上三角矩陣 R
array([[-3.74165739, -4.81070235,  1.87082869],
       [ 0.        , -3.13960871,  7.64426469],
       [ 0.        ,  0.        ,  0.25537696]])
>>> np.dot(Q, Q.T)              # 正交矩陣與其轉置矩陣之內積必為單元矩陣 I
array([[ 1.00000000e+00, -4.66827300e-17,  8.68933848e-17],
       [-4.66827300e-17,  1.00000000e+00,  7.00177831e-17],
       [ 8.68933848e-17,  7.00177831e-17,  1.00000000e+00]])
>>> B=linalg.inv(Q)             # 求正交矩陣 Q 的反矩陣 B
>>> B                                     # 正交矩陣的反矩陣即為本身之轉置
array([[-0.80178373, -0.53452248, -0.26726124],
       [ 0.59152048, -0.77352678, -0.22750788],
       [-0.08512565, -0.34050261,  0.93638218]])
>>> np.allclose(A, np.dot(Q, R))      # 檢查 A 與 QR 是否相等
True  

從三階方陣更能清楚看出正交矩陣的特性. 


6. 用 lu() 函數對方陣進行 LU 分解 :

LU 分解是將一個方陣 A 分解為一個下三角矩陣 L 與上三角矩陣 U 的內積 : 



在線性代數中 LU 分解是高斯消去法的矩陣形式, 也是求反矩陣與計算行列式值的關鍵步驟, 主要用來解線性聯立方程組. 不過, 並非每個方陣都存在 LU 分解, 但若透過一個置換矩陣 (permutation matrix) P 先做行列順序調換, 則任何方陣都可以進行 LU 分解了 : 



所謂置換矩陣是指每列與每行都只有一個 1, 其餘元素均為 0 的方陣, 參考 :

 
SciPy 的 linalg 所提供 lu() 函數的就是做 PLU 分解, 它會傳回 P, L, U 組成之元組 :

P, L, U=linalg.lu(A)

例如 : 

>>> A=np.array([[3, 2, 3], [2, 5, -7], [1, 2, -2]])   # 三階方陣
>>> A   
array([[ 3,  2,  3],
       [ 2,  5, -7],
       [ 1,  2, -2]]) 
>>> P, L, U=linalg.lu(A)        # 對方陣 A 做 PLU 分解
>>> P                                       # 置換矩陣 P 
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])
>>> L                                       # 下三角矩陣 L
array([[1.        , 0.        , 0.        ],
       [0.66666667, 1.        , 0.        ],
       [0.33333333, 0.36363636, 1.        ]])
>>> U                                       # 上三角矩陣 U
array([[ 3.        ,  2.        ,  3.        ],
       [ 0.        ,  3.66666667, -9.        ],
       [ 0.        ,  0.        ,  0.27272727]])
>>> np.allclose(A, P.dot(L.dot(U)))       # 檢查 A 是否等於 PLU 內積
True    

下面範例使用 random 子套件來產生一個隨機三階方陣, 其 LU 分解如下 :

>>> from scipy import random, linalg       # 匯入隨機子套件 random
>>> A=random.randn(3, 3)                          # 產生三階隨機方陣
>>> A     
array([[ 1.05600878,  0.35537349, -0.04829833],
       [ 1.33457398,  0.192862  , -0.03187792],
       [ 0.57623532, -0.36047869, -0.31709815]])
>>> P, L, U=linalg.lu(A)            # 對方陣 A 做 PLU 分解
>>> P                                           # 置換矩陣 P 
array([[0., 0., 1.],
       [1., 0., 0.],
       [0., 1., 0.]])
>>> L                                           # 下三角矩陣 L
array([[ 1.        ,  0.        ,  0.        ],
       [ 0.43177473,  1.        ,  0.        ],
       [ 0.79127032, -0.4569392 ,  1.        ]])
>>> U                                           # 上三角矩陣 U
array([[ 1.33457398,  0.192862  , -0.03187792],
       [ 0.        , -0.44375163, -0.30333407],
       [ 0.        ,  0.        , -0.16167951]])
>>> np.allclose(A, P.dot(L.dot(U)))      
True   

此例的 P 不是單元矩陣, 可見方陣 A 是經過行列順序調換才能做 LU 分解. 

參考 :


2021年10月20日 星期三

Python 學習筆記 : SciPy 的子套件

一直想透過學習 Scipy 把自己虛弱的數學底子好好地補一補, 但時間總是非常有限. 今天得空先將 SciPy 的常用子套件整理成一張表, 常回來看這張表就能提醒自己還有哪些是還沒學的. 

SciPy 系列之前的筆記參考 :


SciPy 是以 Numpy 為基礎建構的開放原始碼科學運算套件, 它依賴底層以 Fortran 實作的演算法可以高效地操作 Numpy 陣列以進行各種數學運算, 其功能足以與商用的 Matlab 或開放原始碼的 Scilab, GSL, 以及 Octave 等軟體相匹敵. 

Scipy 的功能分門別類放在下列子套件中 : 


 SciPy 常用子套件 功能與用途
 scipy.cluster 分群 (clustering) 運算 (例如 K-means)
 scipy.constants 數學與物理常數
 scipy.fftpack 離散傅立葉 (Fourier) 轉換
 scipy.integrate 數值積分與求解常微分方程式
 scipy.interpolate 內插運算 (一維 & 多維內插, 徑向基函數內插, 平滑樣條 spline 等)
 scipy.io 輸入與輸出
 scipy.linalg 線性代數 (求解線性方程組, 特徵值, 特徵向量, 矩陣分解等)
 scipy.signal 訊號處理 (FIR 與 IIR 數位濾波器)
 scipy.stats 統計函數 (連續 & 離散隨機變數的各種機率分布)
 scipy.stats.mstats 遮罩陣列統計函數
 scipy.sparse 稀疏矩陣
 scipy.sparse.linalg 稀疏線性代數
 scipy.sparse.csgraph 壓縮的稀疏圖像函數
 scipy.special 特殊函數
 scipy.spatial 空間演算法與資料結構
 scipy.optimize 最佳化與求根 (求非線性方程組解, 資料擬合, 函數最小值等)
 scipy.odr 正交距離回歸
 scipy.ndimage 多維影像處理 (影像濾波器, 傅立葉轉換, 影像內插, 旋轉等)
 scipy.misc 雜項函數 (求導數, 載入心電圖等)


由於 SciPy 已經相當龐大, 所以後來開發的科學計算相關套件都以獨立形式發布, 而非納入 SciPy 中成為其子套件, 例如機器學習套件 scikit-learn 與 scikit-image 等. 

參考 SciPy 官網教學文件 : 


由於 SciPy 功能收納於各個子套件, 由於有些套件較龐大, 因此很少會用 import scipy 匯入整個套件, 通常會用 from scipy 方式匯入會使用到的子套件 (亦可加上 as 取一個簡名), 例如 : 

from scipy import linalg
from scipy import linalg as la

然後用子套件名稱或簡名呼叫旗下函數來運算, 例如求解線性聯立方程組 :

a_ans=linalg.solve(a_coeff, a_const)     
a_ans=la.solve(a_coeff, a_const)           

雖然可以像 Numpy 那樣直接匯入 SciPy 並取個慣用簡名 sp, 但不建議這麼做 :

import scipy as sp    

但這樣存取函數的路徑就較長了, 例如 : 

a_ans=sp.linalg.solve(a_coeff, a_const)    

另外一種方式是用 from 從子套件直接匯入會用到的函數, 例如 : 

from scipy.linalg import solve

然後直接呼叫該函數做運算 :

a_ans=solve(a_coeff, a_const)

但這種方式要注意命名空間汙染問題 (勿自訂同名函數). 

如果要匯入整個子套件且想要直接呼叫旗下函數, 可在 import 後面使用 *, 例如 : 

from scipy.linalg import *  

科學與工程科系的學生若能好好地掌握 Numpy 與 SciPy 的用法, 在學習研究上必然如虎添翼. 下面是幾本跟 SciPy 相關的好書 :

Python資料運算與分析實戰 (中久喜健司, 旗標, 2018)