解释:def conjugate_gradient(fun, grad, x0, iterations, tol): """ Minimization of scalar function of one or more variables using the conjugate gradient algorithm. Parameters ---------- fun : function Objective function. grad : function Gradient function of objective function. x0 : numpy.array, size=9 Initial value of the parameters to be estimated. iterations : int Maximum iterations of optimization algorithms. tol : float Tolerance of optimization algorithms. Returns ------- xk : numpy.array, size=9 Parameters wstimated by optimization algorithms. fval : float Objective function value at xk. grad_val : float Gradient value of objective function at xk. grad_log : numpy.array The record of gradient of objective function of each iteration. """ fval = None grad_val = None x_log = [] y_log = [] grad_log = [] x0 = asarray(x0).flatten() # iterations = len(x0) * 200 old_fval = fun(x0) gfk = grad(x0) k = 0 xk = x0 # Sets the initial step guess to dx ~ 1 old_old_fval = old_fval + np.linalg.norm(gfk) / 2 pk = -gfk x_log = np.append(x_log, xk.T) y_log = np.append(y_log, fun(xk)) grad_log = np.append(grad_log, np.linalg.norm(xk - x_log[-1:])) gnorm = np.amax(np.abs(gfk)) sigma_3 = 0.01 while (gnorm > tol) and (k < iterations): deltak = np.dot(gfk, gfk) cached_step = [None] def polak_ribiere_powell_step(alpha, gfkp1=None): xkp1 = xk + alpha * pk if gfkp1 is None: gfkp1 = grad(xkp1) yk = gfkp1 - gfk beta_k = max(0, np.dot(yk, gfkp1) / deltak) pkp1 = -gfkp1 + beta_k * pk gnorm = np.amax(np.abs(gfkp1)) return (alpha, xkp1, pkp1, gfkp1, gnorm) def descent_condition(alpha, xkp1, fp1, gfkp1): # Polak-Ribiere+ needs an explicit check of a sufficient # descent condition, which is not guaranteed by strong Wolfe. # # See Gilbert & Nocedal, "Global convergence properties of # conjugate gradient methods for optimization", # SIAM J. Optimization 2, 21 (1992). cached_step[:] = polak_ribiere_powell_step(alpha, gfkp1) alpha, xk, pk, gfk, gnorm = cached_step # Accept step if it leads to convergence. if gnorm <= tol: return True # Accept step if sufficient descent condition applies. return np.dot(pk, gfk) <= -sigma_3 * np.dot(gfk, gfk) try: alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \ _line_search_wolfe12(fun, grad, xk, pk, gfk, old_fval, old_old_fval, c2=0.4, amin=1e-100, amax=1e100, extra_condition=descent_condition) except _LineSearchError: break # Reuse already computed results if possible if alpha_k == cached_step[0]: alpha_k, xk, pk, gfk, gnorm = cached_step else: alpha_k, xk, pk, gfk, gnorm = polak_ribiere_powell_step(alpha_k, gfkp1) k += 1 grad_log = np.append(grad_log, np.linalg.norm(xk - x_log[-1:])) x_log = np.append(x_log, xk.T) y_log = np.append(y_log, fun(xk)) fval = old_fval grad_val = grad_log[-1] return xk, fval, grad_val, x_log, y_log, grad_log
时间: 2024-03-07 19:51:20 浏览: 205
这是一个使用共轭梯度算法来最小化一个标量函数的函数。它的输入包括一个目标函数,一个梯度函数,一个初始值,迭代次数和容差。返回值包括参数的估计值,目标函数在该点的值,目标函数在该点的梯度值以及每次迭代中目标函数梯度的记录。这个算法使用了Polak-Ribiere共轭梯度法,并且在每次迭代中使用了Wolfe条件进行线搜索。算法的大致步骤是:首先计算梯度,然后使用Polak-Ribiere方向来计算下一步方向,然后使用Wolfe条件进行线搜索来确定步长。重复这个过程直到满足停止条件。
相关问题
IFFT_modulation = zeros(symbols_per_carrier + 1, IFFT_bin_length); IFFT_modulation(:,carriers) = complex_carrier_matrix; IFFT_modulation(:,conjugate_carriers) = conj(complex_carrier_matrix);
这段代码的作用是将复数载波矩阵映射到IFFT变换的频域矩阵中。具体来说,假设有symbols_per_carrier个子载波,那么IFFT_bin_length就是整个频域矩阵的长度,包括了零频率分量和负频率分量。carriers是一个长度为symbols_per_carrier的向量,它表示需要映射到IFFT变换频域矩阵中的子载波位置。conjugate_carriers则是对应的共轭位置。
首先,将IFFT_modulation初始化为一个大小为(symbols_per_carrier + 1) x IFFT_bin_length的零矩阵。然后,将复数载波矩阵complex_carrier_matrix映射到IFFT_modulation矩阵的对应子载波位置carriers中,也就是将其赋值给IFFT_modulation矩阵的第carriers个列向量。接着,将复数载波矩阵的共轭值映射到IFFT_modulation矩阵的共轭子载波位置conjugate_carriers中,也就是将其赋值给IFFT_modulation矩阵的第conjugate_carriers个列向量的共轭值。最终,IFFT_modulation矩阵中的所有非零位置都被赋值为复数载波矩阵中对应位置的值和其共轭值的和。这个矩阵将用于进行IFFT变换,生成时域信号。
class SVDRecommender: def __init__(self, k=50, ncv=None, tol=0, which='LM', v0=None, maxiter=None, return_singular_vectors=True, solver='arpack'): self.k = k self.ncv = ncv self.tol = tol self.which = which self.v0 = v0 self.maxiter = maxiter self.return_singular_vectors = return_singular_vectors self.solver = solver def svds(self, A): if self.which == 'LM': largest = True elif self.which == 'SM': largest = False else: raise ValueError("which must be either 'LM' or 'SM'.") if not (isinstance(A, LinearOperator) or isspmatrix(A) or is_pydata_spmatrix(A)): A = np.asarray(A) n, m = A.shape if self.k <= 0 or self.k >= min(n, m): raise ValueError("k must be between 1 and min(A.shape), k=%d" % self.k) if isinstance(A, LinearOperator): if n > m: X_dot = A.matvec X_matmat = A.matmat XH_dot = A.rmatvec XH_mat = A.rmatmat else: X_dot = A.rmatvec X_matmat = A.rmatmat XH_dot = A.matvec XH_mat = A.matmat dtype = getattr(A, 'dtype', None) if dtype is None: dtype = A.dot(np.zeros([m, 1])).dtype else: if n > m: X_dot = X_matmat = A.dot XH_dot = XH_mat = _herm(A).dot else: XH_dot = XH_mat = A.dot X_dot = X_matmat = _herm(A).dot def matvec_XH_X(x): return XH_dot(X_dot(x)) def matmat_XH_X(x): return XH_mat(X_matmat(x)) XH_X = LinearOperator(matvec=matvec_XH_X, dtype=A.dtype, matmat=matmat_XH_X, shape=(min(A.shape), min(A.shape))) #获得隐式定义的格拉米矩阵的低秩近似。 eigvals, eigvec = eigsh(XH_X, k=self.k, tol=self.tol ** 2, maxiter=self.maxiter, ncv=self.ncv, which=self.which, v0=self.v0) #格拉米矩阵有实非负特征值。 eigvals = np.maximum(eigvals.real, 0) #使用来自pinvh的小特征值的复数检测。 t = eigvec.dtype.char.lower() factor = {'f': 1E3, 'd': 1E6} cond = factor[t] * np.finfo(t).eps cutoff = cond * np.max(eigvals) #获得一个指示哪些本征对不是简并微小的掩码, #并为阈值奇异值创建一个重新排序数组。 above_cutoff = (eigvals > cutoff) nlarge = above_cutoff.sum() nsmall = self.k - nlarge slarge = np.sqrt(eigvals[above_cutoff]) s = np.zeros_like(eigvals) s[:nlarge] = slarge if not self.return_singular_vectors: return np.sort(s) if n > m: vlarge = eigvec[:, above_cutoff] ularge = X_matmat(vlarge) / slarge if self.return_singular_vectors != 'vh' else None vhlarge = _herm(vlarge) else: ularge = eigvec[:, above_cutoff] vhlarge = _herm(X_matmat(ularge) / slarge) if self.return_singular_vectors != 'u' else None u = _augmented_orthonormal_cols(ularge, nsmall) if ularge is not None else None vh = _augmented_orthonormal_rows(vhlarge, nsmall) if vhlarge is not None else None indexes_sorted = np.argsort(s) s = s[indexes_sorted] if u is not None: u = u[:, indexes_sorted] if vh is not None: vh = vh[indexes_sorted] return u, s, vh def _augmented_orthonormal_cols(U, n): if U.shape[0] <= n: return U Q, R = np.linalg.qr(U) return Q[:, :n] def _augmented_orthonormal_rows(V, n): if V.shape[1] <= n: return V Q, R = np.linalg.qr(V.T) return Q[:, :n].T def _herm(x): return np.conjugate(x.T) 将上述代码修改为使用LM,迭代器使用arpack
将代码修改为使用LM和迭代器使用arpack的方法如下:
```python
class SVDRecommender:
def __init__(self, k=50, ncv=None, tol=0, which='LM', v0=None, maxiter=None, return_singular_vectors=True, solver='arpack'):
self.k = k
self.ncv = ncv
self.tol = tol
self.which = which
self.v0 = v0
self.maxiter = maxiter
self.return_singular_vectors = return_singular_vectors
self.solver = solver
def svds(self, A):
if self.which == 'LM':
largest = True
elif self.which == 'SM':
largest = False
else:
raise ValueError("which must be either 'LM' or 'SM'.")
if not (isinstance(A, LinearOperator) or isspmatrix(A) or is_pydata_spmatrix(A)):
A = np.asarray(A)
n, m = A.shape
if self.k <= 0 or self.k >= min(n, m):
raise ValueError("k must be between 1 and min(A.shape), k=%d" % self.k)
if isinstance(A, LinearOperator):
if n > m:
X_dot = A.matvec
X_matmat = A.matmat
XH_dot = A.rmatvec
XH_mat = A.rmatmat
else:
X_dot = A.rmatvec
X_matmat = A.rmatmat
XH_dot = A.matvec
XH_mat = A.matmat
dtype = getattr(A, 'dtype', None)
else:
if n > m:
X_dot = X_matmat = A.dot
XH_dot = XH_mat = _herm(A).dot
else:
XH_dot = XH_mat = A.dot
X_dot = X_matmat = _herm(A).dot
dtype = A.dtype
def matvec_XH_X(x):
return XH_dot(X_dot(x))
def matmat_XH_X(x):
return XH_mat(X_matmat(x))
XH_X = LinearOperator(matvec=matvec_XH_X, dtype=dtype, matmat=matmat_XH_X, shape=(min(A.shape), min(A.shape)))
# 获得隐式定义的格拉米矩阵的低秩近似。
eigvals, eigvec = eigsh(XH_X, k=self.k, tol=self.tol ** 2, maxiter=self.maxiter, ncv=self.ncv, which=self.which, v0=self.v0, mode=self.solver)
# 格拉米矩阵有实非负特征值。
eigvals = np.maximum(eigvals.real, 0)
# 使用来自pinvh的小特征值的复数检测。
t = eigvec.dtype.char.lower()
factor = {'f': 1E3, 'd': 1E6}
cond = factor[t] * np.finfo(t).eps
cutoff = cond * np.max(eigvals)
# 获得一个指示哪些本征对不是简并微小的掩码,
# 并为阈值奇异值创建一个重新排序数组。
above_cutoff = (eigvals > cutoff)
nlarge = above_cutoff.sum()
nsmall = self.k - nlarge
slarge = np.sqrt(eigvals[above_cutoff])
s = np.zeros_like(eigvals)
s[:nlarge] = slarge
if not self.return_singular_vectors:
return np.sort(s)
if n > m:
vlarge = eigvec[:, above_cutoff]
ularge = X_matmat(vlarge) / slarge if self.return_singular_vectors != 'vh' else None
vhlarge = _herm(vlarge)
else:
ularge = eigvec[:, above_cutoff]
vhlarge = _herm(X_matmat(ularge) / slarge)
if self.return_singular_vectors != 'u':
u = _augmented_orthonormal_cols(ularge, nsmall) if ularge is not None else None
else:
u = None
if self.return_singular_vectors != 'vh':
vh = _augmented_orthonormal_rows(vhlarge, nsmall) if vhlarge is not None else None
else:
vh = None
indexes_sorted = np.argsort(s)
s = s[indexes_sorted]
if u is not None:
u = u[:, indexes_sorted]
if vh is not None:
vh = vh[index_sorted]
return u, s, vh
def _augmented_orthonormal_cols(U, n):
if U.shape[0] <= n:
return U
Q, R = np.linalg.qr(U)
return Q[:, :n]
def _augmented_orthonormal_rows(V, n):
if V.shape[1] <= n:
return V
Q, R = np.linalg.qr(V.T)
return Q[:, :n].T
def _herm(x):
return np.conjugate(x.T)
```
阅读全文