GNU Octave多数据集全局非线性拟合报错求助（Matlab可正常运行）

问题现象 #

在GNU Octave中执行含共享参数的多数据集全局非线性拟合（部分参数固定）时，出现如下错误：

error: operator *: nonconformant arguments (op1 is 34x1, op2 is 4x1)

相同代码在Matlab中可正常运行。

代码及数据说明 #

拟合代码 #

clear
close all
clc
pkg load optim
D = dlmread('hd', ';'); % 原始数据
bkg = D(1,2:end); % 4个传感器的背景值
x = D(2:end,1); % 输入信号
Y = D(2:end,2:end); % 4个传感器的响应数据
W = 1./Y; % 权重
b0 = [7 .04 .01 .1 .5 2 1]; % 拟合初始猜测值
%% 模型函数
F = @(b) ((bkg + (b(1) - bkg).*(1-exp(-(b(2:5).*x).^b(6))).^b(7)) - Y) .* W;
opts = optimset("Display", "iter");
lb = [5 .001 .001 .001 .001 .01 1];
ub = [];
[b, resnorm, residual, exitflag, output, lambda, Jacob] = ...
lsqnonlin(F,b0,lb,ub,opts)

参数规则 #

参数数组b0中：

• b0(1)、b0(6)、b0(7)为4个数据集的共享参数
• b0(2:5)为各数据集的独有参数

直接在命令行调用F(b0)可得到与Y尺寸一致的34×4矩阵，符合lsqnonlin的输入要求。

原始数据（保存为hd文件，分隔符为;） #

0,0.3105,0.31342,0.31183,0.31117
0.013229,0.329,0.3295,0.332,0.372
0.013229,0.328,0.33,0.33,0.373
0.021324,0.33,0.3305,0.33633,0.399
0.021324,0.325,0.3265,0.333,0.397
0.037763,0.33,0.3255,0.34467,0.461
0.037763,0.327,0.3285,0.347,0.456
0.069405,0.338,0.3265,0.36533,0.587
0.069405,0.3395,0.329,0.36667,0.589
0.12991,0.357,0.3385,0.41333,0.831
0.12991,0.358,0.3385,0.41433,0.837
0.25368,0.393,0.347,0.501,1.302
0.25368,0.3915,0.3515,0.498,1.278
0.51227,0.458,0.3735,0.668,2.098
0.51227,0.47,0.3815,0.68467,2.124
1.0137,0.61,0.4175,1.008,3.357
1.0137,0.599,0.422,1,3.318
2.0162,0.89,0.5335,1.645,5.006
2.0162,0.872,0.5325,1.619,4.938
4.0192,1.411,0.716,2.674,6.595
4.0192,1.418,0.7205,2.691,6.766
8.0315,2.34,1.118,4.195,7.176
8.0315,2.33,1.126,4.161,6.74
16.04,3.759,1.751,5.9,7.174
16.04,3.762,1.748,5.911,7.151
32.102,5.418,2.942,7.164,7.149
32.102,5.406,2.941,7.164,7.175
64.142,7.016,4.478,7.174,7.176
64.142,7.018,4.402,7.175,7.175
128.32,7.176,6.078,7.175,7.176
128.32,7.175,6.107,7.175,7.173
255.72,7.165,7.162,7.165,7.165
255.72,7.165,7.164,7.166,7.166
511.71,7.165,7.165,7.165,7.165
511.71,7.165,7.165,7.166,7.164

解决方法 #

Octave与Matlab的数组广播机制存在差异：b(2:5)是1×4向量，x是34×1向量，Matlab会自动广播为34×4矩阵，但Octave需要显式处理维度匹配问题。可通过以下三种方式修改模型函数：

方法1：使用bsxfun实现元素级乘法 #

F = @(b) ((bkg + (b(1) - bkg).*(1-exp(-(bsxfun(@times, b(2:5), x)).^b(6))).^b(7)) - Y) .* W;

方法2：利用维度转置实现广播 #

F = @(b) ((bkg + (b(1) - bkg).*(1-exp(-(b(2:5).' * x.').^b(6))).^b(7)) - Y) .* W;

方法3：显式扩展维度 #

F = @(b) ((bkg + (b(1) - bkg).*(1-exp(-(repmat(b(2:5), length(x), 1) .* repmat(x, 1, 4))).^b(6))).^b(7)) - Y) .* W;

修改后重新运行代码，即可解决维度不匹配的报错。

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