如何用Householder变换对对称矩阵三对角化？代码通用问题求助

我正在参考一份线性代数笔记，使用Householder变换对对称矩阵进行三对角化，笔记第2页有相关示例。以下是我编写的Python代码：

尽管这段代码对我测试的某个特定哈密顿量有效，但无法通用。我已添加注释，方便理解。据我所知，若矩阵非对称，运行这段代码会得到上三角矩阵。

import numpy as np
from numpy.linalg import matrix_power,norm,eigvals

# 先生成一个维度为d的随机矩阵
d=4
a=np.random.rand(d,d)
# 转为对称矩阵，我研究的是物理中的厄米矩阵
h=(a+a.T)/2
ei,ev=np.linalg.eig(h)
hp=h.copy()

for i in range(d-2):  
    # 记录要置零的元素的范数，存到对角旁的矩阵元素中
    if norm(h[i+2:,i])>0: 
        h[i][i+1]=norm(h[i+2:,i])
       
    for j in range(i+2,d):
        h[i][j]=0
   
   
    w= (hp[i]-h[i])
    
    if norm(w)==0:
        w=w
    else:
        w=w/norm(w)
    p=(np.eye(d) - 2*np.outer(w, w))
   
    h=(p@hp@p.T).round(4)
    
    hp=h.copy()
    

print(h)

内容的提问来源于stack exchange，提问作者sreeram pg