我对python和numpy相对较新，并且正在尝试使用频谱聚类对具有浮点数且尺寸为256x256的密集矩阵进行聚类。由于亲和度矩阵的大小为65536x65536，因此无法计算完整的亲和度矩阵（由于内存限制）。因此，我目前正在计算给定矩阵条目与其5x5局部邻域之间的亲和力，并构建一个稀疏图（以3元组表示）。
为此，我正在使用for循环（基本上是一种滑动的寡妇方法），我认为这并不是最有效的方法。
import numpy as np
def getAffinity(f1, f2):
return np.exp(-np.linalg.norm(np.absolute(f1 - f2))/ 2.1)
G = np.arange(256*256).reshape((256,256))
dim1 = 256 # Dimension 1 of matrix
dim2 = 256 # Dimension 1 of matrix
values = np.zeros(1623076, dtype=np.float32) # To hold affinities
rows = np.zeros(1623076, dtype=np.int32) # To hold row index
cols = np.zeros(1623076, dtype=np.int32) # To hold column index
index = 0 # To hold column index
for i in range(dim1):
for j in range(dim2):
current = G[i, j]
for k in range(np.maximum(0, i-2), np.minimum(dim1 , i+3)): # traverse rows
for l in range(np.maximum(0, j-2), np.minimum(dim2 , j+3)): # traverse columns
rows[index] = i*d1 + j
cols[index] = k*d1 + l
values[index] = getAffinity(current, G[k, l])
index += 1
我想知道是否还有其他有效的方法可以实现相同的目标。
最佳答案
这是一种稀疏矩阵方法。它比循环代码快800倍以上。
import numpy as np
from scipy import sparse
from time import perf_counter as pc
T = []
T.append(pc())
def getAffinity(f1, f2):
return np.exp(-np.linalg.norm(np.absolute(f1 - f2))/ 2.1)
G = 2*np.arange(256*256).reshape((256,256))
dim1 = 256 # Dimension 1 of matrix
dim2 = 256 # Dimension 1 of matrix
values = np.zeros(1623076, dtype=np.float32) # To hold affinities
rows = np.zeros(1623076, dtype=np.int32) # To hold row index
cols = np.zeros(1623076, dtype=np.int32) # To hold column index
index = 0 # To hold column index
for i in range(dim1):
for j in range(dim2):
current = G[i, j]
for k in range(np.maximum(0, i-2), np.minimum(dim1 , i+3)): # traverse rows
for l in range(np.maximum(0, j-2), np.minimum(dim2 , j+3)): # traverse columns
rows[index] = i*dim1 + j
cols[index] = k*dim1 + l
values[index] = getAffinity(current, G[k, l])
index += 1
T.append(pc())
affs_OP = sparse.coo_matrix((values,(rows,cols))).tocsr()
import scipy.sparse as sp
def getAffinity(f1, f2): # similar to @PaulPanzer, I don't think OP is right
return np.exp(-np.abs(f1 - f2)/ 2.1)
def affinity_block(dim = 256, dist = 2):
i = np.arange(-dist, dist+1)
init_block = sp.dia_matrix((np.ones((i.size, dim)), i), (dim, dim))
out = sp.kron(init_block, init_block).tocoo()
out.data = getAffinity(Gf[out.row], Gf[out.col])
return out
T.append(pc())
Gf = G.ravel()
offsets = np.concatenate((np.mgrid[1:3,-2:3].reshape(2,-1).T,np.mgrid[:1,1:3].reshape(2,-1).T), axis=0)
def make_diag(yo,xo):
o = 256*yo+xo
diag = np.exp(-np.abs(Gf[o:]-Gf[:-o])/2.1)
if xo>0:
diag[:xo-256].reshape(-1,256)[:,-xo:] = 0
elif xo<0:
diag[:xo].reshape(-1,256)[:,:-xo] = 0
diag[xo:] = 0
return diag
diags = [make_diag(*o) for o in offsets]
offsets = np.sum(offsets*[256,1], axis=1)
affs_pp = sparse.diags([*diags,[np.ones(256*256)],*diags],np.concatenate([offsets,[0],-offsets]))
T.append(pc())
affs_df = affinity_block()
T.append(pc())
print("OP: {:.3f} s convert OP to sparse matrix: {:.3f} s pp {:.3f} s df: {:.3f} s".format(*np.diff(T)))
diff = affs_pp-affs_OP
diff *= diff.sign()
md = diff.max()
print(f"max deviation pp-OP: {md}")
print(f"number of different entries pp-df: {(affs_pp-affs_df).nnz}")
样品运行：
OP: 23.392 s convert OP to sparse matrix: 0.020 s pp 0.025 s df: 0.093 s
max deviation pp-OP: 2.0616356788405454e-08
number of different entries pp-df: 0
有点解释，一维首先使其保持简单。让我们想象一个实际滑动的窗口，因此我们可以将时间用作直观的轴：
space
+-------------->
|
t | xo... x: window center
i | oxo.. o: window off center
m | .oxo. .: non window
e | ..oxo
| ...ox
v
这里的时间实际上等于空间，因为我们以恒定的速度运动。现在我们可以看到所有窗口点都可以描述为三个对角线。偏移量为0、1和-1，但请注意，由于亲和力是对称的，而0的亲和力很小，因此我们只需要为1计算它们。
现在让我们跳到2D，我们可以做的最小的例子是4x4阵列中的3x3窗口。在行专业这看起来像。
xo..oo..........
oxo.ooo.........
.oxo.ooo........
..ox..oo........
oo..xo..oo......
ooo.oxo.ooo.....
.ooo.oxo.ooo....
..oo..ox..oo....
....oo..xo..oo..
....ooo.oxo.ooo.
.....ooo.oxo.ooo
......oo..ox..oo
........oo..xo..
........ooo.oxo.
.........ooo.oxo
..........oo..ox
相关偏移量是（0,1），（1，-1），（1,0），（1,1）或行主要0x4 + 1 = 1，1x4-1 = 3，1x4 + 0 = 4， 1x4 + 1 =5。另外请注意，这些对角线大多数都不完整，丢失的位由行主要环绕来解释，即在z = y，xx = 3时，右邻z + 1实际上不是右邻y， x + 1;相反，由于行跳动，它是y + 1,0。上面代码中的if-else子句将每个对角线的右边位变为空白。
@DanielF的策略相似，但是利用了图中明显的块结构。
xo.. oo.. .... ....
oxo. ooo. .... ....
.oxo .ooo .... ....
..ox ..oo .... ....
oo.. xo.. oo.. ....
ooo. oxo. ooo. ....
.ooo .oxo .ooo ....
..oo ..ox ..oo ....
.... oo.. xo.. oo..
.... ooo. oxo. ooo.
.... .ooo .oxo .ooo
.... ..oo ..ox ..oo
.... .... oo.. xo..
.... .... ooo. oxo.
.... .... .ooo .oxo
.... .... ..oo ..ox
关于python - 仅计算大型图中局部邻域的亲和力的有效方法，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/57797555/