局部一致概念因子分解（LCCF）是一种强大的无监督聚类算法，它在概念因子分解（CF）的框架下引入了流形正则项，能够在核空间中学习局部一致的低维表示。相比传统NMF，LCCF的基向量是数据点的线性组合（X U），这使得它天然适合核化扩展，能够处理非线性分布的数据，同时通过图正则保持样本间的局部几何结构。
今天分享的这个函数LCCF_Multi正是LCCF算法的核心乘性迭代实现。它基于核矩阵K进行优化，支持归一化拉普拉斯、多次随机初始化选优、灵活的列归一化策略，以及稳定的收敛判断机制，工程实现非常精炼高效。
乘性更新的数学推导
LCCF的目标函数为：
||K - K U V^T||_F² + α Tr(V^T L V)其中K是核矩阵，L是图拉普拉斯矩阵。
通过标准乘性规则推导，可得以下更新公式：
V ← V ⊙ (K U + W V) / (V U^T K U + D V)
U ← U ⊙ (K V) / (K U V^T V)
代码中巧妙地将正则项融入分子和分母，避免显式计算大矩阵乘积，同时利用K的对称性优化计算顺序。
函数实现亮点
高效的乘性更新：
先计算公共小矩阵如U^T K U、V^T V（规模k²，k通常很小）。
正则项以W V和D V形式直接加入分子分母，计算量仅为O(n²k)。
使用元素逐点操作（.* ./）实现更新，保证非负性。
核空间列归一化：
支持对V列L2归一化，或对U列在核空间归一化（sqrt(U^T K U)），确保数值稳定。
每次满足最小迭代后重新归一化。
多次随机初始化：
默认重复10次随机非负初始化，取目标函数最优结果。
有效避免局部最优，提升聚类稳定性。
收敛机制灵活：
使用滑动平均目标函数变化率（meanFit）作为主要停止准则。
支持固定最大迭代次数。
图正则预处理：
支持归一化拉普拉斯（D^{-1/2} L D^{-1/2}），提升正则项效果。
alpha已在前置函数中缩放处理。
目标函数精确计算：
利用trace性质高效计算重构误差：trace(K) - 2 trace(V^T K U) + trace((U^T K U) (V^T V))
分离计算NMF项和拉普拉斯项，便于监控。
使用建议
核选择：高斯核常用于文档/图像数据，线性核退化为线性LCCF。
alpha调参：从10到1000尝试，值越大越强调局部一致性。
k选择：通常设为聚类数或略大。
这个乘性更新模块是整个LCCF算法的性能核心，它将复杂的核化流形学习转化为简单高效的迭代过程。如果你正在实现或复现局部一致概念因子分解，这个函数将提供可靠、高效的优化支持，尤其在文档聚类任务中表现突出。
完整代码（带中文功能注释）
function[U_final,V_final,nIter_final,objhistory_final]=LCCF_Multi(K,k,W,options,U,V)% LCCF_Multi: 局部一致概念因子分解乘性更新规则核心实现%% 输入:% K - 核矩阵 (nSmp x nSmp)，对称正定% k - 隐因子（概念）数量% W - 样本间邻接矩阵% options - 参数结构体（由LCCF传入）% U, V - 初始矩阵（可为空）%% 输出:% U_final, V_final - 最终关联矩阵和概念系数矩阵% nIter_final - 迭代次数% objhistory_final - 目标函数历史differror=options.error;maxIter=options.maxIter;nRepeat=options.nRepeat;minIterOrig=options.minIter;minIter=minIterOrig-1;meanFitRatio=options.meanFitRatio;alpha=options.alpha;Norm=2;NormV=1;nSmp=size(K,1);ifalpha<0alpha=0;endW=alpha*W;DCol=full(sum(W,2));D=spdiags(DCol,0,speye(size(W,1)));L=D-W;% 可选归一化拉普拉斯ifisfield(options,'NormW')&&options.NormW D_mhalf=DCol.^-.5;tmpD_mhalf=repmat(D_mhalf,1,nSmp);L=(tmpD_mhalf.*L).*tmpD_mhalf';clear D_mhalf tmpD_mhalf;L=max(L,L');endifisempty(U)U=abs(rand(nSmp,k));V=abs(rand(nSmp,k));elsenRepeat=1;end[U,V]=NormalizeUV(K,U,V,NormV,Norm);selectInit=1;ifnRepeat==1selectInit=0;minIterOrig=0;minIter=0;ifisempty(maxIter)[obj_NMFhistory,obj_Laphistory]=CalculateObj(K,U,V,L);objhistory=obj_NMFhistory+obj_Laphistory;meanFit=objhistory*10;elseifisfield(options,'Converge')&&options.Converge[obj_NMFhistory,obj_Laphistory]=CalculateObj(K,U,V,L);objhistory=obj_NMFhistory+obj_Laphistory;endendelseifisfield(options,'Converge')&&options.Convergeerror('Not implemented!');endendtryNo=0;whiletryNo<nRepeat tmp_T=cputime;tryNo=tryNo+1;nIter=0;maxErr=1;while(maxErr>differror)% 更新VKU=K*U;% K * UUKU=U'*KU;% U^T * K * UVUKU=V*UKU;% V * U^T * K * Uifalpha>0WV=W*V;DV=repmat(DCol,1,k).*V;KU=KU+WV;% 正则分子项VUKU=VUKU+DV;% 正则分母项endV=V.*(KU./max(VUKU,1e-10));% 乘性更新Vclear WV DV KU UKU VUKU;% 更新UKV=K*V;% K * VVV=V'*V;% V^T * VKUVV=K*U*VV;% K * U * V^T * VU=U.*(KV./max(KUVV,1e-10));% 乘性更新Uclear KV VV KUVV;nIter=nIter+1;ifnIter>minIter[U,V]=NormalizeUV(K,U,V,NormV,Norm);% 重新归一化ifselectInit[obj_NMFhistory,obj_Laphistory]=CalculateObj(K,U,V,L);objhistory=obj_NMFhistory+obj_Laphistory;maxErr=0;elseifisempty(maxIter)[obj_NMF,obj_Lap]=CalculateObj(K,U,V,L);newobj=obj_NMF+obj_Lap;objhistory=[objhistory newobj];meanFit=meanFitRatio*meanFit+(1-meanFitRatio)*newobj;maxErr=(meanFit-newobj)/meanFit;elseifisfield(options,'Converge')&&options.Converge[obj_NMF,obj_Lap]=CalculateObj(K,U,V,L);newobj=obj_NMF+obj_Lap;objhistory=[objhistory newobj];endmaxErr=1;ifnIter>=maxIter maxErr=0;ifisfield(options,'Converge')&&options.Convergeelseobjhistory=0;endendendendendend% 多轮初始化选最优iftryNo==1U_final=U;V_final=V;nIter_final=nIter;objhistory_final=objhistory;elseifobjhistory(end)<objhistory_final(end)U_final=U;V_final=V;nIter_final=nIter;objhistory_final=objhistory;endendifselectInitiftryNo<nRepeat% 重新随机初始化U=abs(rand(nSmp,k));V=abs(rand(nSmp,k));[U,V]=NormalizeUV(K,U,V,NormV,Norm);elsetryNo=tryNo-1;minIter=0;selectInit=0;U=U_final;V=V_final;objhistory=objhistory_final;meanFit=objhistory*10;endendendnIter_final=nIter_final+minIterOrig;Norm=2;NormV=0;[U_final,V_final]=NormalizeUV(K,U_final,V_final,NormV,Norm);%==========================================================================function[obj_NMF,obj_Lap,dV]=CalculateObj(K,U,V,L,deltaVU,dVordU)if~exist('deltaVU','var')deltaVU=0;endif~exist('dVordU','var')dVordU=1;enddV=[];UK=U'*K;% U^T * KUKU=UK*U;% U^T * K * UVUK=V*UK;% V * U^T * KVV=V'*V;% V^T * Vobj_NMF=sum(diag(K))-2*sum(diag(VUK))+sum(sum(UKU.*VV));ifdeltaVUifdVordU dV=V*UKU-UK'+L*V;elsedV=(VUK-K)'*V;endendobj_Lap=sum(sum((L*V).*V));%==========================================================================function[U,V]=NormalizeUV(K,U,V,NormV,Norm)k=size(U,2);ifNorm==2ifNormV norms=max(1e-15,sqrt(sum(V.^2,1)))';V=V*spdiags(norms.^-1,0,k,k);U=U*spdiags(norms,0,k,k);elsenorms=max(1e-15,sqrt(sum(U.*(K*U),1)))';% 核空间L2范数U=U*spdiags(norms.^-1,0,k,k);V=V*spdiags(norms,0,k,k);endelseifNormV norms=max(1e-15,sum(abs(V),1))';V=V*spdiags(norms.^-1,0,k,k);U=U*spdiags(norms,0,k,k);elsenorms=max(1e-15,sum(U.*(K*U),1))';U=U*spdiags(norms.^-1,0,k,k);V=V*spdiags(norms,0,k,k);endendend
