在 Matlab 中分离“纠缠”向量

Question

我有一组三个向量（存储在一个 3xN 矩阵中），它们是“纠缠的”（例如，第二行中的某个值应该在第三行中，反之亦然）。这种“纠缠”是基于查看绘制 alpha2 的图形。为了分离向量，我使用基于差异的方法，计算一个值相对于下一个三个值的差异（例如，将 (1,i) 与 (:,i+1) 进行比较）。然后我取最小值并存储它。该方法可以分离三个向量中的两个，但不是最后一个。

我想知道你们是否可以与我分享您的想法如何解决这个问题（如果可能的话）。我在下面添加了我的代码。

提前致谢！

数字问题：

clear all; close all; clc;
%%
alpha2 = [-23.32    -23.05  -22.24  -20.91  -19.06  -16.70  -13.83  -10.49  -6.70;
    -0.46     -0.33 0.19    2.38    5.44    9.36    14.15   19.80   26.32;
    -1.58     -1.13 0.06    0.70    1.61    2.78    4.23    5.99    8.09];
%%% Original
figure()
hold on
plot(alpha2(1,:))
plot(alpha2(2,:))
plot(alpha2(3,:))

%%% Store start values
store1(1,1) = alpha2(1,1);
store2(1,1) = alpha2(2,1);
store3(1,1) = alpha2(3,1);

for i=1:size(alpha2,2)-1
    for j=1:size(alpha2,1)
        Alpha1(j,i) = abs(store1(1,i)-alpha2(j,i+1));
        Alpha2(j,i) = abs(store2(1,i)-alpha2(j,i+1));
        Alpha3(j,i) = abs(store3(1,i)-alpha2(j,i+1));

        [~, I] = min(Alpha1(:,i));
        store1(1,i+1) = alpha2(I,i+1);

        [~, I] = min(Alpha2(:,i));
        store2(1,i+1) = alpha2(I,i+1);

        [~, I] = min(Alpha3(:,i));
        store3(1,i+1) = alpha2(I,i+1);
    end
end
%%% Plot to see if separation worked
figure()
hold on
plot(store1)
plot(store2)
plot(store3)

Answer 1

通过polyfit 使用外推解决方案：

这个想法非常简单：遍历所有位置 i 并使用 polyfit 将次数为 d 的多项式拟合到从 F(:,i-(d+1)) 到 F(:,i) 的 d+1 值。使用这些多项式推断函数值F(:,i+1)。然后计算最适合这些外推的真实值F(:,i+1) 的排列。如果只涉及几个功能，这应该工作得很好。当然还有一些改进的余地，但对于您的简单设置应该就足够了。

function F = untangle(F, maxExtrapolationDegree)
%//   UNTANGLE(F) untangles the functions F(i,:) via extrapolation.
if nargin<2
    maxExtrapolationDegree = 4;
end
extrapolate = @(f) polyval(polyfit(1:length(f),f,length(f)-1),length(f)+1);
extrapolateAll = @(F) cellfun(extrapolate, num2cell(F,2));
fitCriterion = @(X,Y) norm(X(:)-Y(:),1);

nFuncs = size(F,1);
nPoints = size(F,2);
swaps = perms(1:nFuncs);
errorOfFit = zeros(1,size(swaps,1));
for i = 1:nPoints-1
    nextValues = extrapolateAll(F(:,max(1,i-(maxExtrapolationDegree+1)):i));
    for j = 1:size(swaps,1)
        errorOfFit(j) = fitCriterion(nextValues, F(swaps(j,:),i+1));
    end
    [~,j_bestSwap] = min(errorOfFit);
    F(:,i+1) = F(swaps(j_bestSwap,:),i+1);
end

初步解决方案：（不那么漂亮 - 跳过这部分）

这是一个类似的解决方案，它试图将导数的总和最小化到向量值函数F = @(j) alpha2(:,j) 的某种程度。它通过遍历位置i 并检查i 坐标的所有可能排列来获得函数F(1:i) 的最小半范数。 （我现在实际上想知道是否有任何规范的数学方法来定义半范数，以便我们得到预期的结果......我最初打算使用 H^1 和 H^2 半范数，但它们不太有效...)

function F = untangle(F)
nFuncs = size(F,1);
nPoints = size(F,2);
seminorm = @(x,i) sum(sum(abs(diff(x(:,1:i),1,2)))) + ...
                  sum(sum(abs(diff(x(:,1:i),2,2)))) + ...
                  sum(sum(abs(diff(x(:,1:i),3,2)))) + ...
                  sum(sum(abs(diff(x(:,1:i),4,2))));
doSwap = @(x,swap,i) [x(:,1:i-1), x(swap,i:end)];
swaps = perms(1:nFuncs);
normOfSwap = zeros(1,size(swaps,1));
for i = 2:nPoints
    for j = 1:size(swaps,1)
        normOfSwap(j) = seminorm(doSwap(F,swaps(j,:),i),i);
    end
    [~,j_bestSwap] = min(normOfSwap);
    F = doSwap(F,swaps(j_bestSwap,:),i);
end

用法：

alpha2 = untangle(alpha2); 命令将解开你的函数：

它甚至应该适用于更复杂的数据，比如这些打乱的正弦波：

nPoints = 100;
nFuncs = 5;
t = linspace(0, 2*pi, nPoints);
F = bsxfun(@(a,b) sin(a*b), (1:nFuncs).', t);
for i = 1:nPoints
    F(:,i) = F(randperm(nFuncs),i);
end

备注：我想如果你已经知道你的函数将是二次函数或其他一些特殊形式，RANSAC 对于更多的函数来说会是一个更好的主意。如果函数没有以相同的 x 值间距给出，这也很有用。