使用SciPy实现无平滑最小二乘优化，解决中子谱解谱问题

我正在用邦纳球谱仪做中子谱解谱，这本质是个欠定问题。目前用scipy.optimize的least_squares方法求解，但结果里谱线被“拉平”、峰被展宽，像是隐含的平滑/连续性约束导致结果偏离最优——从问题本质来说，谱线必须存在窄峰才能保证计数准确性。

我的实现代码如下：

import numpy as np
from scipy.optimize import least_squares

def objective(spectrum, kernel, counts, energies):
    residual = np.matmul(kernel, spectrum) - counts
    return residual

def L2_regularization_alogrithm(alpha, beta, kernel, guess, counts, energies):
    guess_L2_reg = guess.reshape((len(guess),))
    counts_L2_reg = counts.reshape((len(counts),))
    used_kernel = kernel.copy()

    kernel_T = np.transpose(used_kernel)
    K = np.dot(kernel_T, used_kernel)
    LM = np.diag(np.diag(K))
    transfer = K + alpha * LM
    counts_L2_reg = np.matmul(kernel_T, counts_L2_reg)

    bounds = (0, np.inf)
    result = least_squares(objective, guess_L2_reg, args=(transfer, counts_L2_reg, energies), bounds=bounds)
    x = result.x
    return x.reshape((len(x), 1))

我试过linear、soft_l1、huber、cauchy和arctan等不同损失函数来缓解问题，当前实现效果最好但仍达不到预期。已知这个问题可解，因为商业软件能输出带清晰窄峰的结果，但无法获取其内部代码。当前最优结果仍存在峰展宽、谱线拉平的问题。

解决思路 #

1. 替换L2正则为L1类稀疏正则
   L2正则天然倾向于生成平滑解，而L1（Lasso）正则能诱导稀疏性，让谱线集中在少数能量点形成窄峰。可以通过增广残差的方式，在least_squares中加入L1正则：

def objective_with_l1(spectrum, kernel, counts, alpha):
    residual = np.matmul(kernel, spectrum) - counts
    # 用增广残差实现L1正则，适配least_squares接口
    return np.concatenate([residual, np.sqrt(alpha) * spectrum])

也可以直接使用scipy.optimize.lbfgsb求解带L1正则的问题，它对稀疏约束的支持更直接。

2. 迭代重加权最小二乘（IRLS）
   针对L1正则，IRLS可将其转化为加权L2问题迭代求解，每轮更新权重放大小谱值的影响，促进稀疏窄峰的形成：

def irls_solve(kernel, counts, guess, alpha, max_iter=10):
    spectrum = guess.copy()
    for _ in range(max_iter):
        # 权重避免除以0，小谱值对应高权重
        weights = 1 / (np.abs(spectrum) + 1e-6)
        weighted_kernel = kernel * weights[np.newaxis, :]
        weighted_counts = counts * weights.flatten()
        # 带非负约束求解加权最小二乘
        result = least_squares(lambda x: np.matmul(weighted_kernel, x) - weighted_counts, 
                              spectrum.flatten(), bounds=(0, np.inf))
        spectrum = result.x.reshape(-1, 1)
    return spectrum

3. 替换正则化方式：二阶差分惩罚
   当前用的alpha*LM是对谱值绝对值的惩罚，换成二阶差分惩罚能精准控制平滑度——它惩罚的是谱线的曲率，而非绝对值，既能抑制无意义噪声，又能保留真实窄峰：

def build_second_diff_matrix(n):
    # 构建二阶差分矩阵，D@spectrum输出相邻点的二阶差
    D = np.zeros((n-2, n))
    for i in range(n-2):
        D[i, i] = 1
        D[i, i+1] = -2
        D[i, i+2] = 1
    return D

def objective_with_smooth(spectrum, kernel, counts, alpha):
    residual = np.matmul(kernel, spectrum) - counts
    smooth_residual = alpha * np.matmul(build_second_diff_matrix(len(spectrum)), spectrum)
    return np.concatenate([residual, smooth_residual])

4. 优化初始猜测
   如果初始猜测是平坦谱，算法容易陷入平滑解的局部最优。可以基于已知的中子峰能量位置，在初始猜测中对应能量点设置较高值，其余设为小值，引导算法向有峰的方向收敛。

5. 使用专业欠定问题求解器
   用cvxpy构建带稀疏约束的优化问题，能更灵活定义目标函数，求解器会自动处理稀疏性和非负约束：

import cvxpy as cp

def cvxpy_sparse_solve(kernel, counts, alpha):
    n = kernel.shape[1]
    spectrum = cp.Variable(n, nonneg=True)
    # 最小化拟合误差+L1稀疏正则
    objective = cp.Minimize(cp.norm(kernel @ spectrum - counts, 2) + alpha * cp.norm(spectrum, 1))
    prob = cp.Problem(objective)
    prob.solve()
    return spectrum.value.reshape(-1, 1)

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