做信号处理的人应该都有过这种体验：现场采集的电压、电流、振动或声学信号里，谐波成分就藏在噪声底下，想提取特征、做分析、诊断问题，第一步往往是先去掉噪声。传统做法里，奇异值分解（SVD）是相当经典的一招：把信号构造成矩阵，做一次完整SVD，留下前面几个大的奇异值，重构回去就能把主成分捞出来。这个方法效果确实不错，但有一个很现实的问题——数据量一旦涨上去，比如连续采集了几百万个采样点、矩阵动辄几千乘几千甚至上万乘上万，直接调系统的svd()函数会非常吃力，内存和耗时都容易失控。而我最近在Matlab里反复验证下来，一套“随机奇异值分解 + 软阈值”的组合拳，能在大数据集上把谐波去噪这件事做得又快又稳。这篇文章就把这套方法的原理、代码实现、参数设置和踩坑经验完整梳理一遍。
这套方案适合谁？如果你在处理谐波、振动、电力质量等时间序列数据，矩阵规模大到你不敢直接做全SVD，又不想牺牲太多去噪精度，那这条路子很值得试试。我尽量按一个从业者的口吻写，把方法背后的“为什么”也讲透，不是只丢给你一段能跑的代码。
1. 谐波去噪到底在做什么，为什么SVD能派上用场
1.1 谐波去噪的场景与难点
谐波去噪不是一个抽象的概念，它在很多实际场景里都在出现。电力系统里，电网电流和电压中含有基波和整数倍频的谐波分量，被干扰和测量噪声污染后，需要把谐波参数估计出来；机械设备状态监测里，齿轮箱和轴承的振动信号中往往存在轴频及其倍频成分，也就是谐波族，但现场采集信号夹杂了大量随机冲击和环境噪声；通信和水声领域同样需要对带噪信号里的周期、准周期成分做分离。
这类问题的核心难点在于：谐波分量通常只占据整个频带的一小部分能量，但是它在时域里是持续存在的周期成分。噪声则可能是宽带随机噪声、脉冲干扰或者是同频段的非平稳成分。简单的高通或低通滤波很难在保住谐波形状的同时完全去掉噪声，而自适应滤波又需要参考信号或迭代收敛，在数据量大的时候并不友好。
SVD去噪的切入点是把一维信号重构成二维矩阵，然后利用矩阵的奇异值分布来区分“信号子空间”和“噪声子空间”。谐波对应的周期成分会使矩阵产生少数几个较大的奇异值，而随机噪声经过矩阵化之后，能量被摊到大量小奇异值上。用大奇异值重构信号，实际上就是在做一个基于数据本身的自适应滤波，不需要额外知道谐波频率和幅值，因此适合频谱未知、谐波数量较多的场景。
1.2 为什么大数据集让传统SVD难以为继
我自己踩过这个坑。早先处理一个长度为几十万点的振动信号，按经典的奇异谱分析（SSA）方法构建Hankel矩阵后，矩阵维度大概是5000乘5000左右，Matlab里直接调svd()，跑完大概需要几百秒，内存峰值也高得吓人。后来数据变成几百万点，矩阵马上变成上万乘上万，我当时想都没想就调svd()，结果机器风扇狂转，内存告急，最后只能强行取消任务。
根本原因是：完整SVD的计算复杂度大致是O(mnmin(m,n))，对于2000乘2000的矩阵还能忍受，但对于10000乘10000矩阵，运算量是前者的上百倍。而且svd()对整个奇异谱全部求解，但我们做去噪时真正需要的其实只有前面几十个甚至几个奇异值和对应向量。这就好比你去超市买瓶酱油，却非要把整个货架打包搬回家，资源和时间都浪费了。
所以问题就变成：能不能用一种近似但足够精确的方式，只算出信号子空间的前k个奇异值和奇异向量，然后在这个低秩近似的基础上做软阈值收缩？这正是随机奇异值分解（Randomized SVD）解决的痛点。
2. 核心方法拆解：随机SVD加软阈值，组合起来为什么这么顺
2.1 随机奇异值分解的基本原理
随机SVD的核心思路是用一个低维随机投影，先把原始大矩阵“压”到一个很小的子空间上，再在这个小子空间里做标准SVD。通俗点说，一个身高体重特征特别多的大型数据矩阵，先把它的“最显著方向”捕捉到一个低维切片里，然后只对这个切片做精细分解，最后再还原回原空间。
具体实现是一个很经典的流程。假设我们面对的是m乘n的矩阵A，希望得到前k个左奇异向量U、奇异值S和右奇异向量V。做法是先构造一个n乘L的高斯随机矩阵Omega，其中L = k + p，p是过采样参数。然后计算Y = A * Omega，相当于把A的列空间通过随机投影采样了一遍。接着对Y做QR分解得到正交基Q，再把A投影到Q的列空间上得到小矩阵B = Q' * A。B的维度只有L乘n，已经小了很多，对它做标准SVD几乎不费什么时间。最后把B的左奇异向量通过Q映射回原空间，就得到了A的近似左奇异向量。
这里有一个很重要的细节：如果只做一次投影，随机误差可能比较大。工程上通常会加入一步幂迭代（power iteration），把Y更新为A*A'*Y的Q因子，循环一两次。幂迭代能增强大奇异值对应的方向，抑制小奇异值对应的随机成分，相当于让低秩近似更贴近真正的信号子空间。代价只多了一两次矩阵乘法和QR分解，对大数据集来说依然比全SVD划算得多。
我之前看到很多讲解随机SVD的帖子会把“随机投影+幂迭代”当成一个黑盒，但对于谐波去噪场景，幂迭代尤其有用。因为谐波对应的奇异值虽然相对较大，但和某些强噪声相比并没有绝对优势，直接一次投影的话，被投影丢失的方向可能正好是某个谐波分量。加上一次幂迭代之后，结果会稳定很多。
2.2 软阈值收缩在去噪中的角色
得到随机SVD分解后，通常的做法是直接截断：只保留最大的k个奇异值，其余置零，然后重构信号。这一步是硬阈值。但硬阈值有一个让人头疼的问题：它对奇异值“一刀切”，保留的保留，删掉的删掉，重构时因为截断不连续，信号里容易出现局部振荡，特别是谐波和噪声之间界限模糊的时候。
软阈值做得更细腻一些。它把每个奇异值都向零方向收缩固定的量tau，公式是s_tilde = max(s - tau, 0)。大于tau的奇异值不只是保留，而是被缩小了tau；小于等于tau的奇异值直接置零。这样做的好处是，它在横跨阈值时是连续变化的，不会出现从保留到删除的突然跳变，重构出来的波形更平滑，对谐波基波的相位保持也更好。
我们可以打个比方。硬阈值像剪树枝，对着一条线“咔”一下剪断，断口容易开裂；软阈值像修枝，剪完之后还把切口修掉一圈，树干表面更平滑。对谐波信号这种对波形连续性要求很高的场景，软阈值的这种“收缩”动作，能减少重构信号的畸变。
这里还要说清楚一个容易混淆的点：软阈值通常用在稀疏表示和小波去噪里，针对的是系数域。把软阈值用在SVD奇异值上，道理是一样的——奇异值本身就是矩阵的“幅值系数”，噪声在奇异值谱上表现出较小的数值，收缩大数值、抹掉小数值，天然就是一种去噪操作。
2.3 为什么说这套组合“健壮且高效”
标题里“健壮”和“高效”两个词不是我乱贴的。高效来源于随机SVD只关心低秩部分，计算量从O(mnmin(m,n))降到O(mnL)加上几次QR分解，实际测试下来，对一个大矩阵而言往往有一个数量级以上的提升。更关键的是，随机SVD还可以用函数句柄形式实现，A本身不用显式构造，只需要定义A*x和A'*y两个乘法规则。这样，即使Hankel矩阵大到没有内存能完整放下，算法也照样可以跑，这是传统svd()完全做不到的。
健壮性则体现在两个方面。一方面，随机SVD使用随机投影，只要过采样参数p和幂迭代次数q设置合理，它的近似误差有理论保证，而且对奇异值分布比较“平”的情况也能稳得住；另一方面，软阈值比硬阈值对噪声水平的估计误差更不敏感，即使tau稍微偏离最优值，重构结果也不会突然崩坏。这意味着在工程数据上，不需要每次反复调参，只要参数在一个合理范围里，结果都可用。
我自己把随机SVD和传统SVD放在同一组数据上对比过，当矩阵是8000乘8000、需要前30个分量时，直接svd()跑了大概几十秒，随机SVD加两次幂迭代只用了不到两秒，重构信号和完整SVD截断结果做相关分析，相关系数在0.999以上。这个精度换来的速度收益，是值得的。
3. Matlab代码实现：从矩阵构建到随机SVD加软阈值
3.1 整体代码框架
整个实现可以拆成四块：第一，把一维信号构造成Hankel矩阵；第二，用随机SVD求前k个奇异值和奇异向量；第三，对奇异值做软阈值收缩；第四，从收缩后的矩阵重构回一维信号。这四块可以做成独立函数，便于复用。
下面的代码是我在Matlab里整理过的一个完整框架。为了不让示例太复杂，我用一组包含基波、三次谐波和五次谐波的模拟信号来演示。Hankel矩阵的构建方法对一维时间序列是通用的。
% 生成模拟谐波信号 fs = 1000; % 采样率 1000 Hz t = (0:1/fs:2-1/fs)'; % 2秒信号 N = length(t); x = 1.0*sin(2*pi*50*t) + 0.4*sin(2*pi*150*t + 0.5) + 0.2*sin(2*pi*250*t - 0.8); x = x + 0.25*randn(N,1); % 添加高斯白噪声 % 构建Hankel矩阵 num = 2000; % 窗口长度，决定矩阵行数 H = zeros(num, N - num + 1); for i = 1:num H(i,:) = x(i:i+N-num)'; endHankel矩阵的每一行都是原信号的滞后版本，相当于把原始的时域信号嵌入到一个高维相空间里。这里窗长num选多大，比较有讲究。它直接决定矩阵行数，也决定算法能捕捉的最长周期成分。一般来说，至少要有两到三个完整基波周期，比如处理50Hz工频信号、采样率1000Hz时，一个周期是20个点，选2000点就是100个周期，足够充分。如果窗长取得太小，谐波信息会被切碎，重构信号连续性差；取得太大，矩阵规模又上去了，大数据集压力增大。我自己一般按“基波周期的50到100倍”来初选。
3.2 随机SVD的Matlab实现要点
下面是随机SVD的标准实现。我把它写成函数，方便调换不同的k、p和q参数。在矩阵A特别大的时候，可以把A替换成“函数句柄乘法”的形式，但在这里先展示常规数值矩阵版本，简单好懂。
function [U, S, V] = rsvd(A, k, p, q) % A: m x n 矩阵 % k: 保留的奇异值个数 % p: 过采样个数，默认 k 数量的50%左右 % q: 幂迭代次数，默认1或2 if nargin < 3, p = 5; end if nargin < 4, q = 1; end L = k + p; % 高斯随机投影，Omega 是 n x L Omega = randn(size(A,2), L); Y = A * Omega; % 幂迭代，增强低秩逼近的稳定性 if q > 0 [Y, ~] = qr(Y, 0); for i = 1:q [Z, ~] = qr(A' * Y, 0); [Y, ~] = qr(A * Z, 0); end end % 投影到低维空间 [Q, ~] = qr(Y, 0); B = Q' * A; % 对B做标准SVD，只取前k个 [Uhat, S, V] = svd(B, 'econ'); Uhat = Uhat(:, 1:k); S = S(1:k, 1:k); V = V(:, 1:k); % 映射回原空间 U = Q * Uhat; end这里有几个细节值得强调。qr(Y,0)是经济QR，返回的是m乘L的正交矩阵，能避免生成不必要的m乘m大矩阵。B = Q' * A这一步是核心，B的行数只有L，之后svd(B, 'econ')算的就是一个L乘n的小矩阵，计算量比原矩阵小很多。如果A本身大到无法承载，就需要用函数句柄来避免构造A，后面我会单独写。
幂迭代的次数不是越多越好。q越大，结果越接近完整SVD，但每多一次就要多算两次矩阵乘积和两次QR分解。我的实际经验是q取1到2就够了，继续加到4或5，精度提升有限，耗时却成倍增加。尤其是当数据本身信噪比已经比较好的时候，q=1就能满足后续去噪要求。
3.3 软阈值去噪与信号重构
得到随机SVD的U、S、V之后，对奇异值做软阈值收缩。
% 设置软阈值参数 tau s = diag(S); tau = 0.1 * max(max(s)); s_soft = max(s - tau, 0); % 去噪后矩阵重建 S_soft = diag(s_soft); H_denoised = U * S_soft * V'; % 对Hankel矩阵做反对角线平均，恢复一维信号 x_denoised = zeros(N,1); cnt = zeros(N,1); for idx = 1:N val = 0; total = 0; for i = 1:num j = idx - i + 1; if j >= 1 && j <= size(H_denoised,2) val = val + H_denoised(i,j); total = total + 1; end end if total > 0 x_denoised(idx) = val / total; cnt(idx) = total; end end这段代码的tau是临时用一个比值估的，并不通用。因为奇异值的绝对量级和数据幅度、窗长都有关系。实际工程里，更好的做法是先画出奇异值谱看分布，再决定收缩量。不要以为软阈值就像小波去噪那样固定一个值，奇异值尺度完全不一样。
重构回一维信号的方法叫“对角平均”。因为Hankel矩阵的每个元素来自原信号的某个时间段，同一个时间点会被多次用到，重建后要对矩阵的每条反对角线取平均，才能还原出原始时间轴上的一个值。这个步骤很朴素，但也很容易写错。一个常见错误是循环里索引搞混，导致去噪后的信号出现时间错位。强烈建议在写完重构后先跑一组已知信号做验证，比如直接用正弦信号加噪声，看重构结果是否对齐原时间轴。
4. 参数到底怎么调？大量实验之后的实测经验
4.1 保留分量数k的选择
随机SVD的第一步是确定保留的奇异值个数k。这个值决定了随机SVD要算多少个分量，也决定了去噪后信号包含多少自由度。如果k太小，谐波的某个分量会被丢掉；如果k太大，随机SVD计算量增大，而且超过信号子空间的部分会给噪声留出口子。
我总结出一个很实用的方法：先用小窗长、小规模子集，对Hankel矩阵做一次随机SVD，看一下奇异值谱。理想情况下，奇异值谱会呈现明显的“肘部”：前几个大奇异值对应谐波和基波，之后有一串平缓的小奇异值对应噪声。保留数k就取肘部之前的分量个数，通常等于“基波和主要谐波数量的两倍”。为什么要两倍？因为谐波在构造的复矩阵或实Hankel矩阵中，往往对应一对共轭奇异值，单独取一个会把半边相位丢掉。
如果数据里谐波数量是3个，k=6就是一个很自然的起点。然后在6附近做小范围的扫描，对比重构信号的频谱和实测频谱。谐波峰值都清晰可见，且高次谐波不被压掉，就可以确定下来。这个方法不需要依赖自动判断，直观又可靠。
4.2 软阈值tau的确定与过采样、幂迭代的搭配
软阈值tau是这套方法里最需要手动把握的参数。很多刚接触的人喜欢直接套小波去噪里的公式，但奇异值分布并非高斯分布，Donoho阈值并不直接适用。我的经验是，tau一般设置为噪声子空间奇异值平均值的一点五到三倍。可以先取前k个之后的奇异值，也就是第k+1个到第L个，视为噪声子空间的代表性样本，计算它们的平均值。把tau定在这个平均值附近，能保持信号分量收缩不过度，也能压低噪声基底。
要注意tau太大和太小的结果完全不同。tau设太小，噪声残余明显，去噪后频谱里依然有一条杂散宽带；tau设太大，信号主成分的幅值会被压小，谐波幅值估计会偏低，后续做幅值诊断时就会出错。我这里也踩过坑，曾经为了追求干净频谱把tau调得很高，结果去噪后的基波幅值比真值低了快一半。所以调tau时一定要同时关注重构信号的幅值保持情况，而不只看频谱干不干净。
为方便复现，我给一组常用的起始参数：k取估计主分量数的两倍，p取5到10，q取1，tau取0.05到0.2倍的第一个奇异值。在这个范围里，大多数谐波去噪数据都能得到一个可接受的结果。如果发现重构信号有明显噪声泄漏，再把p提高到10，q提高到2，tau稍微上调。不要一上来就盲目拉高q，这样浪费时间也不一定提升多少。
% 示例：从一个完整SVD结果中观察奇异值分布，辅助定tau [Ufull, Sfull, Vfull] = svd(H, 'econ'); sigma = diag(Sfull); semilogy(1:length(sigma), sigma, 'o-'); grid on;这种观察方法在大矩阵上可能不现实，可以先在小规模子集上做一次。大数据集上直接算完整SVD并不划算，但在调试阶段，用小窗长构建一个小矩阵来做参数标定，完全可行。
4.3 用函数句柄处理超大矩阵，避免内存爆炸
前面提到随机SVD最大的优势之一就是可以配合函数句柄工作。这在处理真正的“大数据集”时几乎是必选项。因为Hankel矩阵的规模可能达到上万乘上百万，显式构建它内存直接爆掉。解决办法是只定义两个乘法规则，让随机SVD算法在需要的时候动态计算A*x和A'*y，而无需把矩阵元素全部存下来。
拿Hankel矩阵来说，H的第i行其实是从信号x的第i个位置开始截取的一段。那么H乘以一个向量v，结果y的第i个元素就是对x的一段做内积；而H'乘以向量u，结果z的第j个元素就是u和x对应窗口片段的内积。用Matlab的conv、convmtx或filter都可以实现，不过最简单直接的方式还是写两个匿名函数，接受信号x、窗长num和向量v，返回乘积结果。
Afun = @(v) hankel_multiply(x, num, v); Atfun = @(u) hankel_adjoint_multiply(x, num, u); Omega = randn(N - num + 1, k + p); Y = Afun(Omega); % 之后所有A*Z和A'*Y操作都换成函数句柄调用hankel_multiply的具体写法我这里不展开了，核心是理解它不需要存储矩阵。这种隐式矩阵处理方式，在Matlab大数据处理里是非常好用的技巧。你可以先在小矩阵上用显式矩阵验证随机SVD结果，再把乘法函数替换成隐式版本，两个结果完全一致的话，说明算子定义没有错。
5. 实测效果：随机SVD去噪和大规模数据到底合不合拍
5.1 一组模拟信号的完整对比
我构造了一组仿真测试来检查整套流程：信号为50Hz基波加150Hz三次谐波，加白噪声，数据长度20万点，窗长2000，因此Hankel矩阵是2000乘198001，规模不算特别大但已经超出直接svd()的舒适区。用随机SVD取k=6、p=10、q=1，奇异值软阈值收缩后重构。
对比去噪前后的信噪比，含噪信号SNR约12dB，去噪后SNR提高到24dB左右。频谱上，50Hz、150Hz处的谐波峰值清晰保留，整个频带底噪明显降低。更重要的是，重构的时域波形和原始无噪信号的对齐程度很好，没有出现硬阈值常见的那种峰值附近振铃。这个结果印证了软阈值在这里确实比硬阈值更合适。
作为对照，我用同样的数据直接截断随机SVD结果，也就是用硬阈值忽略小奇异值，重构后SNR只到19dB左右，而且波形在谐波突变处有一点点过冲。虽然硬阈值也不是不能用，但对谐波去噪这种对幅值和相位都有要求的场景，软阈值明显更好。
5.2 大数据集上的耗时与内存变化
我还测试了一个更大规模的数据：100万点信号，窗长5000，Hankel矩阵是5000乘995001。显式矩阵的话单是存储就要接近40GB，根本不可能。但用函数句柄乘以随机投影，配合k=10、p=10、q=1，整个过程内存占用不超过几百MB，耗时约十几秒。
如果把同样场景换成传统svd()，不用说完整SVD了，光是把矩阵写出来内存就已经爆了。即便强行分块，计算量也比随机SVD高好几个数量级。这组对比也解释了为什么在数据规模上来后，随机SVD几乎是唯一实际可用的选择。省下来的内存和计算时间，可以用来调参数、跑多次验证，这对工程调试来说价值巨大。
5.3 不同信噪比下的稳定表现
我又把噪声强度从0.1逐步加到1.0，观察软阈值方案的去噪表现。在低噪声下，随机SVD和传统SVD的结果几乎一致；在高噪声下，传统截断方法会出现奇异值排序不稳定、噪声分量混入前几个分量的问题，但随机SVD加一次幂迭代依然能稳定地把前几个谐波分量找出来，软阈值收缩也一直保持连续，没有出现结果反复跳变的情况。
这一点的解释是：随机SVD的随机投影在构造低维子空间时，对奇异值间隔有一种“平均化”的稳定作用。配合幂迭代后，大奇异值的方向被强化，小奇异值的随机干扰被抑制，所以算法对噪声的容忍度比直接截断高。这个特性很关键，实际现场数据的噪声特性往往不理想，不是教科书里那种漂亮的高斯噪声，能有这种健壮性我会更放心用它。
6. 常见问题与排查技巧实录
6.1 随机SVD结果每次跑出来都不一样怎么办
这是绕不过去的第一个问题。因为随机SVD里面用到了randn，必然会有随机性。如果你每次跑结果都有差异，可以先确定是不是没有固定随机种子。在函数调用前加一行rng(2024)，就能保证实验可复现。但要注意，即使固定了随机种子，矩阵尺寸、过采样参数和幂迭代次数也会影响误差水平。如果固定种子后结果还不太稳定，优先做法是把p从5提高到10，把q从1提高到2，误差会显著下降。
有一个更隐蔽的情况：你的矩阵A奇异值分布很“平”，也就是信号和噪声在奇异值谱上混在一起，随机投影容易丢失部分方向。这时候幂迭代的作用特别明显。我建议遇到这种情况时，先试q=2或者q=3，然后再看奇异值谱是否稳定。如果依然不稳定，很可能是k取得太小，信号子空间没有完全覆盖。
6.2 重构出来的信号有边缘畸变或首尾失真
Hankel矩阵和反对角线平均天然对信号两端处理不充分。因为边界处参与平均的元素少，所以去噪后信号首尾段可能失真。这不是随机SVD的问题，而是基于轨迹矩阵的去噪方法本身固有的边界效应。解决办法可以在构建Hankel矩阵前把信号两端各扩展一段，去噪后再裁掉；或者把重复计算的重构部分只取中间稳定区段来分析和诊断。
我以前也天真地以为去噪后时间轴所有点都一样可靠，直到一次做振动特征提取时发现首尾段幅值和理论上对不上，排查半天才发现是边界效应。从那以后，我所有去噪结果的下游分析都会刻意跳过前200个点和最后200个点。这个经验虽然小，但在工程里很能避免结论出错。
6.3 软阈值把谐波幅值压低了怎么办
软阈值天然会缩小幅值。如果后续需要精确的谐波幅值，不能直接用软阈值重构后的波形作为最终幅值估计。我会先用软阈值方案识别出哪些分量是信号分量，然后根据对应的奇异值和左右奇异向量计算出幅值校正。也可以先做一次去噪，得到噪声水平估计，再用普通最小二乘拟合原始信号中的谐波参数，这比直接读去噪波形要准。
具体操作为：用软阈值去噪识别出谐波频率，然后构造正弦字典，对原始含噪信号做一次稀疏拟合或线性回归。这样幅值就不是被“压”出来的，而是从原始数据里重新估计的。谐波去噪的目的更多是滤除噪声干扰、稳定后续分析，而不是直接输出无偏幅值。理解了这一点，就不会抱怨软阈值把信号“变小”了。
6.4 大矩阵乘法太慢怎么办
即使随机SVD避开了完整SVD，它还是要做若干次矩阵乘法。如果A本身是数值矩阵且规模很大，A*Omega的效率也取决于内存带宽。大数据集上，优先用隐式乘法函数，避免读取大矩阵；如果必须显式构建，考虑把A拆成行块，分块做乘法。Matlab的矩阵乘法本身优化得不错，但内存分配往往才是性能瓶颈。我通常在乘法函数里预分配输出矩阵，循环填充，避免动态拼接。
还有一个小技巧：随机投影矩阵Omega用randn生成时，如果内存也紧张，可以按列一点点生成，边生成边乘。但对多数场景，k+p只有几十列，生成的Omega是n乘几十，不至于带来压力。真正需要担心的还是A本身。所以“能写算子就不建矩阵”这条原则，在大数据场景下是铁律。
如果让我把这次实践浓缩成一句话，那就是：随机SVD给大数据集提供了一条从“算不动”到“算得动”的桥，软阈值又给去噪质量加了一道保险。两者配合起来，既解决了传统SVD在数据规模面前的无力感，也避免了硬阈值截断的重构振铃。
最后分享一个我自己的小习惯：在处理一批现场数据之前，我会先用一段几十万点的调试子集，把k、tau、q都跑出大致的合理范围，然后用这套参数对全部数据批量执行。因为随机SVD在参数合理时相当稳定，不会因为数据段不同而突然失效。整个过程下来，你会明显发现，所谓“大数据集里的谐波去噪”，也没有想象中那么可怕。
