1 目标函数(总)
论文笔记:Temporal Regularized Matrix Factorization forHigh-dimensional Time Series Prediction_UQI-LIUWJ的博客-CSDN博客
1.1 求解W
我们留下含有W的部分:
?
然后对wi求导
线性代数笔记:标量、向量、矩阵求导_UQI-LIUWJ的博客-CSDN博客
而是一个标量,所以放在xi的左边和右边没有影响
所以
也即:
?对应的代码如下:(假设sparse_mat表示 观测矩阵)
from numpy.linalg import inv as inv
for i in range(dim1):
#W矩阵的每一行分别计算
pos0 = np.where(sparse_mat[i, :] != 0)
#[num_obs] 表示i对应的有示数的数量
Xt = X[pos0[0], :]
#[num_obs,rank
vec0 = sparse_mat[i, pos0[0]] @ Xt
#sparse_mat[i, pos0[0]] 是一维向量,
#所以sparse_mat[i, pos0[0]] @ Xt 和 sparse_mat[i, pos0[0]].T @ Xt 是一个意思,
#输出的都是一个一维向量
#[rank,1]
mat0 = inv(Xt.T @ Xt + np.eye(rank))
#[rank,rank]
W[i, :] = mat0 @ vec0
?其中:
| vec0 = sparse_mat[i, pos0[0]] @ Xt | | mat0 = inv(Xt.T @ Xt + np.eye(rank)) |
1.2 求解X
我们留下含有X的部分
表示逐元素乘积 (两个向量a和b,ab可以用diag(a) b表示)
当t=1~ld的时候,我们没有什么事情,所以此时我们更新X的方式和之前的W差不多
同理,X的更新方式为:
而当t≥ld+1的时候,我们就需要考虑了
对于任意xt(我们令其为xo),他会出现在哪些中呢?
首先 是?
对xo求导,有:
其次,是所有的?
对每一个l,有用的项就是xo相关的项,于是我们可以写成,对每一个l的
对xo求导,有
还要注意一点,此时我们考虑的是下标为o+l的AR,所以o+l需要大于ld,小于T,也就是此时o的范围是o>ld-l ,<T-l【也就是说,在ld之前的xo也会有一定的AR的更新项】
于是我们可以写成
几部分拼起来,有
=0
=
+
所以xo(o≥ld+1)的更新公式为
[+]
?代码如下:
for t in range(2):
pos0 = np.where(sparse_mat[:, t] != 0)
#[num_obs] 表示i对应的有示数的数量
Wt = W[pos0[0], :]
#[num_obs,rank],t可见的wi
Mt = np.zeros((rank, rank))
#\sum diad (theta * theta)的那个矩阵
Nt = np.zeros(rank)
#Nt 相当于是 sum theta* x_{t-l}
if t < np.max(time_lags):
#这一个if,else考虑到的是首先的部分
Pt = np.zeros((rank, rank))
#t<ld的时候,是没有λx I的
Qt = np.zeros(rank)
#t<ld的时候,没有过去时间段的回归项
else:
Pt = np.eye(rank)
#t>ld+1的时候 有λx I
Qt = np.einsum('ij, ij -> j', theta, X[t - time_lags, :])
#theta [d,rank]
#X[t - time_lags, :] [d,rank]
'''
对于每一个theta和X[t - time_lags, :]中的j (j∈range(rank))
我们输出是一个rank长的向量,第j维是每个theta的(i,j)元素
和对应的X的(i,j)元素相乘,然后求和
每个theta的(i,j)元素表示向量第j个元素的第i个time lag的AR权重
每个X的(i,j)元素表示向量第j个元素的第i个time lag的数据
相乘再1求和,就是输出的第j个元素,有之前的时序数据加权求和得到的自归回值
'''
'''
换一个角度理解,就是theta和X[t - time_lags, :]逐元素乘积,
再从一个[d,rank]矩阵,压缩成一个rank长的向量
即np.sum(theta*X[t - time_lags, :],axis=0)
'''
if t < dim2 - np.min(time_lags):
#这个if 考虑的是其次的部分,所以需要t+l>ld 且t+l<T
#(ld是max(time_lags),T是dim2
#t+min(lag)小于T
if t >= np.max(time_lags) and t < dim2 - np.max(time_lags):
index = list(range(0, d))
#t>=ld,同时t+max(ld)也没有超过最大的时间限制,此时所有的d个time lad都可以使用
else:
index = list(np.where((t + time_lags >= np.max(time_lags))
& (t + time_lags < dim2)))[0]
#在两头,计算可以用的o+l
for k in index:
Ak = theta[k, :]
#[rank]
Mt += np.diag(Ak ** 2)
#对应的是Σdiag(θk*θk)
theta0 = theta.copy()
theta0[k, :] = 0
#第k行赋值为0的作用,就是Σ_{l'∈L-{l}}
Nt += np.multiply(Ak,
X[t + time_lags[k], :] - np.einsum(
'ij, ij -> j',
theta0,
X[t + time_lags[k] - time_lags, :]))
'''
AK——θl [rank]
X[t + time_lags[k], :] x_{o+l}
np.einsum('ij, ij -> j',theta0, X[t + time_lags[k] - time_lags, :])
对于每一个theta0和XX[t + time_lags[k] - time_lags, :]中的j (j∈range(rank))
我们输出是一个rank长的向量,第j维是每个theta的(i,j)元素
和对应的X的(i,j)元素相乘,然后求和
每个theta的(i,j)元素表示向量第j个元素的第i个time lag的AR权重
每个X的(i,j)元素表示向量第j个元素的第i个time lag的数据
相乘再求和,就是输出的第j个元素,有之前的时序数据加权求和
得到的自归回值(之前的时间序列不包括l,因为那一行被赋值为0了)
'''
#也就是Nt相当于其次括号里面的第一项
vec0 = Wt.T @ sparse_mat[pos0[0], t] + lambda_x * Nt + lambda_x * Qt
#Wt.T @ sparse_mat[pos0[0], t] [rank]
#lambda_x * Nt 第一项
#lambda_x * Qt 第二项
mat0 = inv(Wt.T @ Wt + lambda_x * Mt + lambda_x * Pt + lambda_x * eta * np.eye(rank))
#分别是第一项、第四项、第二项、第三项
X[t, :] = mat0 @ vec0
?vec0:
?+]
mat0:
3 更新θ
我们留下和θ (θk)有关的部分
关于θk求导
4 总结
x:
t ∈ 1~ld:
t ≥ld+1???[+]
?
|