2021年10月31日 星期日

2021 年第 44 周記事

本周回頭讀 Python 進階課題例如迭代器, 產生器, 修飾器等, 以及 zip, enumerate 等內建函數, 我之前對此不求甚解, 一知之半解湊合著用, 但這樣在後續學習中經常會踢鐵板, 所以回頭顧根本. 最近在想我好像輕重緩急搞錯了, 機器學習對我而言雖然最有吸引力, 但似乎不是現在最迫切的, 應該先做對的事, 不要操之過急. 

本周跟水某看完一部 2018 年的舊片 "雖然 30 但仍 17", 由藝術專業戶申惠善主演, 她演技真的很好, 動作與表情都很貼切. 其實我比較喜歡專業管家 Jennifer, 連水某都說好想有個 Jennifer (先中樂透再說啦), 哈哈. 朴恩彬演的戀慕本周已到第七集, 快到一半了 (20 集), 原來者隱君早就知道世子是女的. 

本來不想種菜了, 但今天去種子行買兩包南瓜子, 又順帶買了十幾株甜玉米與三株絲瓜, 所以下午又掄起鋤頭在原來種地瓜葉的地方整理出一畦地來種玉米, 晚上兩個手肘果然酸死了, 畢竟已經半年沒下田了. 下午小舅來菜園整理小番茄, 說絲瓜適合夏天種, 現在不太好種, 對耶, 好像買錯了, 應該買大白菜才對. 

2021年10月30日 星期六

市圖還書 3 本 (工程數學 & 邏輯設計)

這三本書是上學期幫二哥借的, 課已修完先還回去了 : 
我有想要復習工程數學 (for 數位信號處理), 但不是現在, 先記下來唄. 

媽媽的聲音

今天在 Line TV 看 "請回答 1988" 第十集後半段時, 看到正煥的爸爸金老闆在生日那天都會莫名的憂鬱, 他老婆趁他走開時把放在衣櫥裡不肯丟的雜物拿出去外面的垃圾桶丟掉, 但在其中發現了一捲紀錄正峰與正煥兄弟 9 歲與 3 歲時唱歌的錄音帶, 於是就留下來帶回去大家一起聽, 老婆聽到小孩以前可愛的歌聲不禁笑了起來, 但金老闆卻若有所思, 說我聽見媽媽的聲音了, 那時應該是媽媽幫他們帶小孩吧! 錄音那天剛好是金老闆生日, 媽媽要他唱首歌來聽, 說生你的時候不知有多辛苦 ... 金老闆聽到一半就站起來走到門口, 對著戶外滂沱大雨唱起那首歌 ... 令人好感傷. 


我覺得我跟金老闆個性有點像, 都好念舊 (有紀念價值的東西都不肯丟, 哈哈), 又喜歡搞笑 (但只有成德善會跟金老闆玩這種冷遊戲). 他以前是金星電子的員工, 後來好像被裁員, 在正峰中奧運彩券之前金家很窮, 日子不好過. 金星電子就是今天的 LG (1995 年改名, Lucky Goldstar 的縮寫), 我家第一台錄放影機就是 Goldstar 的, 是以前在電腦公司第一年尾牙抽中的, 真的很 Lucky ㄟ (但最 Lucky 的是第二年老總特別抽出的三兩黃金紀念幣, 超 Lucky, 正面是公司圖樣, 背面刻老闆名字, 一直都給母親保存直到她過世前幾年因怕鄉下小偷而交給我, 金價飆到最高點時媽曾問說要不要賣掉, 我怎麼捨得啊). 

自從數位照相機問世後就很少在洗相片了, 所有的照片與影像都放在硬碟裡, 現在則大都是存在手機中, 但這是有風險的, 因為硬碟可能會壞掉, 換手機也可能無意中忘了保存照片, 還是要找時間挑照片洗出來, 相簿才是最好的時光紀錄器, 但影音檔就只能放硬碟了. 我也是好久沒聽到媽媽的聲音了, 我隔一兩年都會買顆硬碟將以前的照片與錄影錄音檔備份轉存到新的硬碟, 這是有錢也買不到的珍貴資產, 想到時就會點出來看一看. 

2021年10月29日 星期五

買 BRUNO 電子壓力鍋 BOE058-BR

由於福利社點數還剩下 1000 多點快到期了, 看來看去沒甚麼特別想買, 就挑了下面這款 BRUNO 的電子壓力鍋, 價格是 2590 元, 比露天便宜, 扣掉 1155 點數, 需自付 1435 元 : 






這款體積較小, 只有 1.5 公升, 相關介紹影片與調理方法如下 : 









這款是日本品牌但是在強國製造, 品質當然比不上純日本製. 本來想買去年買的鍋寶電子壓力過 (容量較大), 但今年那款沒上架, 反正放在高雄主要是給菁菁做便當菜用, 小一點也無所謂. 

2021年10月27日 星期三

Python 學習筆記 : 可迭代物件 (iterable) 與迭代器物件 (iterator)

在 Python 中最常用的迴圈應該是 for in 語法了, 它是透過拜訪可迭代物件來走訪元素, 特別是 for i in range(n) 這樣的語句最常見, 其中的 range() 函數就是一種可迭代物件 (iterable), 它與迭代器 (iterator) 是 Python 迭代協定 (iteration protocol) 裡面所定義的兩種物件類別 : 
  • Iterable (可迭代物件) :
    實作了可傳回一個迭代器物件的 __iter__() 方法. 
  • Iterator (迭代器物件) : 
    繼承自 Iterable 並實作了 __iter__() 與 __next__() 方法的資料流物件, 前者會傳回迭代器 (iterator) 物件本身; 後者則會傳回資料流中的下一個元素 (若已無下一個元素將觸發 StopIteration 例外).
簡言之, Iterator 是 Iterable 的子類別, 故迭代器物件一定是可迭代物件, 但可迭代物件不一定是迭代器物件. 參考 :


Python 中的容器型別例如字串, 元組, 串列, 與陣列 (有序), 以及字典與集合 (無序) 等資料型態都是可迭代物件 (不是迭代器物件), 它們均繼承自 iterable 與 container 類別, 如下圖所示 : 




與迭代有關的 Python 常用內建函數有下列兩個 : 
  • iter(iterable) : 將傳入的可迭代物件轉成迭代器物件 (添加 __next__() 方法). 
  • next(iterator) : 傳回迭代器物件的下一個元素. 
注意, 由於 Python 的容器物件都不是迭代器, 而是可迭代物件, 因此不可以直接把他們傳入 next() 函數, 應該先用 iter() 轉成迭代器才可以傳給 next(), 例如 :

>>> next('abc')      # 字串是可迭代物件, 不是迭代器物件, 不可直接傳入 next()
Traceback (most recent call last):
  File "<pyshell>", line 1, in <module>
TypeError: 'str' object is not an iterator    
>>> str_itr=iter('abc')     
>>> type(str_itr)    
<class 'str_iterator'>
>>> next(str_itr)    
'a'
>>> next(str_itr)   
'b'
>>> next(str_itr)  
'c'

可見 str 物件傳給 iter() 後會傳回一個 str_iterator 類別的迭代器物件, 這個才能傳給 next(). 連續呼叫 next() 會將迭代器物件中的內容逐一傳回, 這與 for 迴圈得作用一樣. 但是在 for 迴圈中卻不需要用 iter() 先把可迭代物件轉成迭代器, 因為 for 迴圈會自動轉換, 例如 : 

>>> for i in 'abc':      # for 迴圈會自動將 'abc' 轉成迭代器 
    print(i)   
    
a
b
c

當然若要手動轉成迭代器也可以, 但那是多此一舉 :

>>> for i in iter('abc'):    
    print(i)   
    
a
b
c

檢驗一個物件是否為 Iterable 或 Iterator 物件可利用 Python 內建模組 collections 裡的兩個基礎類別 IterableIterator, 它們是 Python 所有容器物件的根類別, 需先用 import 指令匯入 :  

from collections import Iterable, Iterator

然後呼叫內建函數 isinstance() 並將這兩個類別作為第二參數傳入, 根據傳回值是 True 或 False 即可檢驗一個物件是否為一個迭代器或是一個可迭代物件 :

isinstance(my_obj, Iterable)       # 傳回 True 表示 my_obj 為可迭代物件
isinstance(my_obj, Iterator)       # 傳回 True 表示 my_obj 為迭代器物件

另外一個方法是呼叫 hasattr() 函數檢查物件是否有 '__iter__' 與 '__next__'  這兩個字串, 若只有 '__iter__'  就是 iterable 物件, 兩個都有就是 iterator 物件 : 

hasattr(my_obj, '__iter__') and not hasattr(my_obj, '__next__')    # True=可迭代物件
hasattr(my_obj, '__iter__') and hasattr(my_obj, '__next__')           # True=迭代器物件

其實迭代器只要檢查有無 __next__() 方法是否存在即可. 

例如 : 

>>> from collections import Iterable, Iterator      # 匯入迭代根類別
>>> isinstance('abc', Iterable)  
True
>>> isinstance('abc', Iterator)   
False    
>>> hasattr('abc', '__iter__')       # 字串是可迭代物件
True
>>> hasattr('abc', '__next__')      # 字串不是迭代器物件
False    

可見字串只是可迭代物件, 不是迭代器物件. 除了字串以外, 其它四個 Python 容器 (tuple, list, dict, set) 以及 range 物件都是如此.


1. range 物件 : 

呼叫 range() 函數會傳回一個可迭代的 range 物件, 但它不是迭代器, 例如 : 

>>> from collections.abc import Iterable, Iterator 
>>> type(range(5))                               # range() 傳回一個 range 物件
<class 'range'>    
>>> isinstance(range(5), Iterable)      # range 是一個可迭代物件
True
>>> isinstance(range(5), Iterator)      # range 不是一個迭代器
False
>>> for i in range(5):                           # for 迴圈會自動將可迭代物件轉成迭代器
    print(i)   
    
0
1
2
3
4
>>> next(range(5))                               # 必須是迭代器物件才能傳給 next()
Traceback (most recent call last):
  File "<pyshell>", line 1, in <module>
TypeError: 'range' object is not an iterator   
>>> n=iter(range(5))                            # 呼叫 iter() 將 range 物件轉成迭代器
>>> next(n)                                           # 呼叫 next() 會傳回迭代器的下一個元素
0
>>> next(n)   
1
>>> next(n)     
2
>>> next(n)  
3
>>> next(n)
4
>>> next(n)   
Traceback (most recent call last):
  File "<pyshell>", line 1, in <module>
StopIteration
 

2. tuple 物件 : 

元組 (tuple) 物件也是可迭代物件, 不是迭代器物件 :

>>> from collections.abc import Iterable, Iterator
>>> t=(1, 2, 3)     
>>> type(t)     
<class 'tuple'>   
>>> isinstance(t, Iterable)      
True
>>> isinstance(t, Iterator)      
False  
>>> hasattr(t, '__iter__')       
True
>>> hasattr(t, '__next__')       
False   


3. list 物件 : 

串列 (list) 物件也是可迭代物件, 不是迭代器物件 :

>>> from collections.abc import Iterable, Iterator
>>> l=[1, 2, 3]     
>>> type(l)     
<class 'list'>   
>>> isinstance(l, Iterable)      
True
>>> isinstance(l, Iterator)      
False  
>>> hasattr(l, '__iter__')       
True
>>> hasattr(l, '__next__')       
False   


4. dict 物件 : 

字典 (dict) 物件也是可迭代物件, 不是迭代器物件 :

>>> from collections.abc import Iterable, Iterator
>>> d={1: 'a', 2: 'b', 3: 'c'}     
>>> type(d)     
<class 'dict'>   
>>> isinstance(d, Iterable)      
True
>>> isinstance(d, Iterator)      
False  
>>> hasattr(d, '__iter__')       
True
>>> hasattr(d, '__next__')       
False   


5. set 物件 : 

集合 (set) 物件也是可迭代物件, 不是迭代器物件 :

>>> from collections.abc import Iterable, Iterator
>>> s={1, 2, 3}     
>>> type(s)     
<class 'set'>   
>>> isinstance(s, Iterable)      
True
>>> isinstance(s, Iterator)      
False  
>>> hasattr(s, '__iter__')       
True
>>> hasattr(s, '__next__')       
False   

參考 :


Chrome 分頁書籤紀錄 (2021)

因為 Chrome 越來越占用資源, 而我又常在幾個學習項目之間切換, 所以常常開了十幾個分頁, 導至筆電效能降低, 以前都用協作平台紀錄這些分頁等待有空時再看, 現在新的協作平台不好用, 所以改記在此處 : 


兩本 Python 進階的好書 : Python Tricks 與 Python for Geeks

Python 語法簡單易學, 通常初學者一周內就可學會基本用法開始撰寫程式; 但若要寫出簡潔專業的程式碼則還有許多密技需要學習, 下面這兩本書就是邁向 Python 進階的好書. 

1. Python 神乎其技 (2020 全新超譯版) : 


Source : 博客來


此書原文版是知名 Python 線上學習網站 realpython.com 站長 Dan Badar 於 2017 年自行出版的 "Python Tricks", 國內知名電腦書作家江良志於 2018 年將其翻成中文版 (旗標), 後來又增添內容再版稱為超譯版 (但原文並未出新版) : 

Python Tricks: A Buffet of Awesome Python Features (2017)


Source : 天瓏


超譯版添加了 Python 3.6 以後新增的功能, 內容來自 Dan Bader 在其教學網站 realpython.com 的更新. 由於這本書實在寫得太棒了, 有一位資料科學家 "好豪" 將超譯版章節與教學網頁之對應整理在其部落格, 非常值得參考 : 



2. Python for Geeks (Packt 2021) : 


Source : Packt


此書作者雖然已是 Nokia 的技術主管, 但仍不停地寫程式. 並將其多年程式經驗的精華濃縮於這本書, 若能熟練書中所介紹的技巧 (特別是第六章 "Advanced Tips and Tricks in Python"), 讀者將具備開發複雜與大型程式專案的能力. 書中範例程式碼可在 GitHub 下載, 參考 : 


2021年10月25日 星期一

二哥的機率與統計課本

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



Source : 天瓏


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

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

2021年10月24日 星期日

用 Python 監測光世代網路速度 (一)

月初把鄉下的光世代升速後測試網速確實有提高, 平均大約是 40M 左右 : 





這比我的 Samsung Note 8 在鄉下的網速要兩倍, 下面是今天在鄉下測試 4G 行動網路的速度 : 





雖然手機相對比較慢, 但比較穩, 我昨晚回到家下載 1GB 檔案, 早上起來發現下載網路連線中斷, 重新下載不到半小時又中斷, 改用行動網路就 OK, 三小時即下載完畢. 可見光世代雖然是固定網路並不是非常穩, 反而比行動網路空中介面還差. 

我想寫個爬蟲放在鄉下家的樹莓派 Pi 3 上, 固定從下面這個測速網站擷取網速後記錄下來, 然後每天寄一份統計圖給我, 藉以觀察鄉下家光世代的網速狀況, 查詢時間間隔 60~180 秒, 每 5 分鐘取平均網速存入檔案或資料庫 : 


觀察其網頁原始碼, 主要發出 request 的部分為 form2 這張表單 : 





<form name='form2' method='post' action='index_noflash.php'>
  <input type='hidden' name='count_time' value=2>
  <input type='hidden' name='count_sum' value=49>
  <div align="center">
    <table align="center" width='600' align="cener">
      <tr>
        <td width="50%" align="right">請選擇要測的檔案大小:</td>
        <td></td>
      </tr>
      <tr>
         <td align="right">
           <select name='size'>
             <option value='4' selected='selected'>4MB</option>
             <option value='8'> 8MB</option>
             <option value='16'>16MB</option>
             <option value='24'>24MB</option>
             <option value='32'>32MB</option>
             <option value='56'>56MB</option>
             <option value='64'>64MB</option>
             <option value='128'>128MB</option>
             <option value='256'>256MB</option>
           </select>
         </td>
         <td align="left">
            <input type="submit" name="Submit2" value=" 再次測試 ">
         </td>
       <t/r>
     </table>
   </div>
</form>

按下 "再次測試" 時會送出三個參數 : 
  • count_time
  • count_sum
  • size




資料先記錄下來, 有空再來測試看看. 

2021-10-25 補充 :

今天搜尋 PyPi 發現 speedtest-cli 這個套件 :

https://pypi.org/project/speedtest-cli/ 

看來只要用這個套件就能完成任務, 應該不需要寫爬蟲了. 另外下面這個 tespeed 也可以, 但似乎是 Python 2 的, 年代有點久遠 :

測試網速的Python工具
https://github.com/Janhouse/tespeed

料理實驗 : 三杯杏鮑菇

今天去市場買菜時老闆說香菜太貴要用買的, 就塞了一堆九層塔給我 (其實我菜園有很多), 想說丟掉又可惜, 所以中午就拿來煮三唄杏鮑菇, 參考食譜如下 : 


材料 :
  • 杏鮑菇 3~5 朵
  • 九層塔 適量
  • 薑 一段
  • 蔥 兩條
  • 蒜 4~5 片
  • 辣椒 適量
  • 糖 適量
  • 醬油 適量
作法 : 
  1. 杏鮑菇洗淨切塊備用. 蔥切段, 薑切片, 蒜片壓碎備用.
  2. 起油鍋, 放入蒜, 薑片, 蔥白炒香, 放入杏鮑菇, 加入少許鹽, 拌炒後蓋上鍋蓋悶兩分鐘, 杏鮑菇會出水. 
  3. 加入兩茶匙糖與適量醬油拌炒, 待水將收乾, 加入米酒, 九層塔, 少許麻油, 拌炒至汁收乾即可起鍋. 



嗯, 還蠻下飯的. 

2021 年第 43 周記事

好快, 一周又過去了, 節氣已來到霜降. 週四去打 AZ 第二劑, 總算打好打滿, 診所醫師說下個月可以打流感疫苗, 雖然說好, 但我一向對打病毒類疫苗興趣缺缺, 藥也是一樣, 能不吃就不吃. 

本周學習進度有點慢, 大概是追劇追太兇了, 哈. 剛看完韓素希演的 "以吾之名", 只有七集, 男主居然 GG, 讓人傻眼. 這是韓版無間道, 黑白互相臥底, 整部戲除了結尾外我覺得拍得還不錯, 韓素希為這部戲應該吃足苦頭, 戲中她倒吊拉單槓是真人嗎? 為戲能練到這樣真是令人佩服啊! 最近中午吃飯時間看的是 "請回答 1988", 這齣每集都很長, 吃完飯只看到一半, 要兩天才看完一集, 總集數 20 集, 雖然很想知道結局 (已猜到德善嫁給正煥, 善宇則姊弟戀娶了寶拉), 但我打算慢慢看. 周二晚上是看朴恩彬與路雲演的 "戀慕", 難得的古裝劇. 

鄉下家的小雖貓的兒子雖然還在吃奶, 但看來已經會跟媽媽搶吃乾糧了, 紙箱已關不住它, 只好讓它在車庫自由行動, 找不到時只要問小雖貓妳兒子咧? 它就會去把小喵喵找出來, 通常都躲在車庫牆邊的鐵架下面. 昨晚看她們母子倆在曬穀場邊打鬥嬉鬧還真有趣, 我想小雖貓應該是同時在教她兒子怎麼狩獵吧. 

今天去市場發現菜價又貴了一些, 香菜與蔥都要買的, 而且不便宜. 打算本周菜園整一整, 下周開始來種下列作物 : 
  • 絲瓜
  • 茄子
  • 香菜
  • 甜玉米
但菜園長滿冬馬齊, 要清除這種可惡的雜草恐怕要花不少時間. 

2021年10月22日 星期五

打 AZ 疫苗第二劑

昨天 10/21 早上請公假去黃昭文診所打第二劑 AZ (這次預約時沒看到巨蛋選項, 全都是診所), 距離 7/29 打第一劑已過 11 周, 上回看到一篇新聞說 AZ 間隔長一些效果較好, 本想下一輪再打, 但因為已收到可以預約的簡訊, 那就打吧! 因為愛心那邊訪視也要求打完兩劑. 

上回打 AZ 第一劑除了第二天上班累累想睡覺外, 並無發燒頭痛等副作用, 這次為了避免下午昏昏欲睡下午也請了假, 但整個下午精神卻挺好, 居然把 SciPy 的線性代數部分都大致測完, 哈. 

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)