MATLAB的答案

使曲线拟合得更快

110(30天)
Marek Dostal
Marek Dostal 2018年11月24日
编辑: 布鲁诺陈德良 2022年11月10日
你好,
我想让我的曲线拟合利用GPU更快。是可能的吗?如果是如何?
谢谢
1评论
马特·J
马特·J 2019年3月25日
Marek问题描述复制以下供参考:
我的大脑核磁共振数据(每个主题都有大约。0、5毫升压)与10个不同的参数,我想在每个立体像素符合这一10分由特定biexponencial曲线(所以我大脑接收参数的地图)。之前我试图优化脚本使用并行计算(大脑1 ~ 4小时),但我只用parfor(~ 1, 25小时),因为我还是初学者Matlab我很好奇如果其他用户有经验GPU,能帮我把这个比parfor快。
谢谢
Marek

登录置评。

在回答这个问题。

接受的答案

马特·J
马特·J 2018年11月26日
编辑:马特·J 2019年6月21日
我想在每个立体像素符合这一10分由特定biexponencial曲线(所以我大脑接收参数的地图)
@Marek,
不符合一个立体像素在一个循环中。相反,编写模型作为一个N维曲线其中N是体素和应用的数量 lsqcurvefit 这个模型。不要让 lsqcurvefit 使用其默认的有限差分计算导数。相反,提供自己的雅可比矩阵x *操作符使用 JacobianMultiplyFcn选项 ( 编辑: 或者提供一个稀疏的块对角雅可比矩阵计算)。
你的问题就是一个很好的例子 我讨论与偶像 :模型函数和 JacobianMultiplyFcn 可能是GPU-optimized如果 lsqcurvefit 可以使用gpuArray变量没有经常做GPU / CPU转移。然而,我预计仍将相当快如果你使用适当的Matlab向量化技术。
16条评论
显示 15年长的评论
莎拉·戈德堡
莎拉·戈德堡 2019年5月6日
嗨,马特,
我一直试图概括和雅可比矩阵的解决方案,与雅可比矩阵计算/计算/不仅雅可比矩阵乘法计算。雅可比矩阵相乘方法是目前不工作。
这里k 1 d符合n个点的数量和p参数,n = n x k, p = x k。雅可比矩阵大小确实是(n x p)——我检查lsqcurvefit的输出。
我遇到以下问题与雅可比矩阵乘法函数(myJ):输入向量的大小Y(资料片)或(Px1) 或(Px2)。 这怎么可能?
下面的代码。我只是看
大小(Y)
谢谢,
莎拉
清晰的所有;关闭所有;文件关闭所有;clc;
s = rng (“洗牌”);
rng (828016308,“旋风”);
全球历史;
历史= [];
%的阴谋
doPlot = false;%真实;
num_plot = 4;%多少符合情节
%合适参数
alg =“trust-region-reflective”;%’trust-region-reflective’或‘levenberg-marquardt”
displaySetting =“关闭”;%的iter”、“了”、“最后”,等等。
麦克斯特= 1 e5;
maxFunEvals = 1 e5;
tolfun = 1 e-6;
%的选项
选择= optimoptions (“lsqcurvefit”,“麦克斯特”麦克斯特,“MaxFunEvals”maxFunEvals,“OutputFcn”@myoutput,“显示”displaySetting,“算法”alg,“tolfun”,tolfun);
optsJ = optimoptions (“lsqcurvefit”,“麦克斯特”麦克斯特,“MaxFunEvals”maxFunEvals,“OutputFcn”@myoutput,“显示”displaySetting,“算法”alg,“tolfun”tolfun,的雅可比矩阵,“上”);
optsJm = optimoptions (“lsqcurvefit”,“麦克斯特”麦克斯特,“MaxFunEvals”maxFunEvals,“OutputFcn”@myoutput,“显示”displaySetting,“算法”alg,“tolfun”tolfun,“SpecifyObjectiveGradient”,真的,“JacobianMultiplyFcn”@(动力系统,Y,国旗,xdata) myJ(动力系统,Y,标志));
%的模型。生成函数处理评估f (f_func)及其偏导数(dfdp_func)。
[f_func, dfdp_func num_params] = sym_model_derivatives ();
%一般2 d模型
num_g = 2;%的高斯函数
npoints = (3、5);%图像大小[r、c] = [y、x], ROI故意不对称进行调试
min_width = 1;%所以我们没有σ= = 0
a0 = 30;
a_std = 1;
mu_std = 0.1;
sigma0 = (2, 1);%小心,不要让这些太小了!
sigma_std = 0.1;
θ= 0 *(π/ 180);%的弧度
theta_std = 10 *(π/ 180);
bg = 10;
bg_std = 5;
mu0 = npoints / 2;%高斯函数大致在中心
我= 0 (npoints);
%生成数据:num_g的图像大小(im)
ims = 0 (num_g、大小(im, 1),大小(im, 2));%图像组大小
亩= mu0 ' * 1 (1, num_g) + mu_std * randn (2, num_g);% 1 /尺寸
(= sigma0 ' * 1 (1, num_g) + sigma_std * randn (2, num_g);% 1 /尺寸
当= a0 + a_std * randn (1, num_g);% 1 /高斯
(((< 0.1)= min_width;%,以避免除0
θ=θ+ theta_std * randn (1, num_g);
英国地质调查局= bg + bg_std * randn (1, num_g);
(X, Y) = meshgrid(1:尺寸(im, 2), 1:尺寸(im, 1));% x - y坐标——奇怪的秩序,这样可以重塑正确的数据
x_ind = X (:);% x,列向量
y_ind = Y (:);% y,列向量
为毫米= 1:num_g
%参数顺序:%幅度%行(y) %坳(x) %θ3 * 3 * % Sigma_row (y) % Sigma_col (x) %的背景
c0 (mm:) =((毫米)亩(2毫米)亩(1毫米)θ(毫米)(开关(开关)(2毫米)(1毫米)英国地质调查局(毫米)];
ims (mm:) = model_func (c0 (mm,:), x_ind, y_ind, f_func);每2 d图% 1维向量
结束
%的猜测
猜测= 0 (num_g num_params);
%参数顺序:%幅度%行% %θ3 * 3 * % Sigma_row上校% Sigma_col %的背景
为毫米= 1:num_g
[maxval ~] = max (ims (mm,:));
tmp_im =重塑(ims (mm,:),大小(im, 1),大小(im, 2));
[xpos, ypos] =找到(tmp_im = = maxval);
[~,half_x] = min (abs (tmp_im (ypos:) - maxval / 2));
[~,half_y] = min (abs (tmp_im (:, xpos) - maxval / 2));
猜测(mm,:) = (maxval xpos ypos 0 (abs (half_x (1) -xpos) + min_width) (abs (half_y (1) ypos) + min_width) 0];%添加的东西很小,以避免0
结束
% % N x 2 d -设置
im_data = 0 (num_g prod(大小(im)));
xdata = kron ((1, num_g), x_ind) ';
ydata = kron ((1, num_g), y_ind) ';
guesses_ND = 0 (num_g num_params);
%的图片
为毫米= 1:num_g
im_data (mm,:) = ims (mm,:);
结束
%的猜测
guesses_ND =猜测;
% % N x 2 d
历史= [];
results_ND = 0 (num_g num_params);
流(1,' \ n ']);
抽搐
xdata f_ND = @ (c) model_func_ND (c xdata ydata f_func);%需要通过x_ind, y_ind model_func
[results_ND, resnormN residualN、exitflagN outputN, lambdaN, jacobianN] = lsqcurvefit (f_ND、guesses_ND xdata, im_data,[],[],选择);%每一行是1对象
toc
流(1,“Nx2D \ n”]);
resnormN
historyN =历史;
exitflagN
% % N x 2 d +雅可比矩阵
results_ND_J = 0 (num_g num_params);
流(1,' \ n ']);
抽搐
f_ND = @ (c, xdata) model_func_ND_J (c xdata ydata、f_func dfdp_func);%需要传递参数x_ind, model_func y_ind等等
[results_ND_J, resnormN_J residualN_J、exitflagN_J outputN_J, lambdaN_J, jacobianN_J] = lsqcurvefit (f_ND、guesses_ND xdata, im_data, [], [], optsJ);%每一行是1对象
toc
流(1,“Nx2D +雅可比矩阵\ n”]);
resnormN_J
exitflagN_J
% % N x 2 d +雅可比矩阵相乘
历史= [];
results_ND_J = 0 (num_g num_params);
流(1,' \ n ']);
抽搐
f_ND = @ (c, xdata) model_func_ND_Jm (c xdata ydata、f_func dfdp_func);%需要传递参数x_ind, model_func y_ind等等
[results_ND_Jm, resnormN_Jm residualN_Jm、exitflagN_Jm outputN_Jm, lambdaN_Jm, jacobianN_Jm] = lsqcurvefit (f_ND、guesses_ND xdata, im_data, [], [], optsJm);%每一行是1对象
toc
流(1,“Nx2D +雅可比矩阵相乘\ n”]);
resnormN_Jm
historyNm =历史;
exitflagN_Jm
返回
% % - - - - - - - - - - - - -功能
%模型:
% f =振幅*。* exp (- (A * (x-Col)。^ 2 + 2 * B * (x-Col)。* (y-Row) + C * (y-Row)。^ 2))+背景
% = cos(θ)^ 2 / (2 * Sigma_col) + sin(θ)^ 2 / (2 * Sigma_row)
% B =罪(2 *θ)* (1 / Sigma_row ^ 2 - 1 / Sigma_col ^ 2) / 4;
% C =罪(θ)^ 2 / (2 * Sigma_col) + cos(θ)^ 2 / (2 * Sigma_row);
%参数顺序:%幅度%行% %θ3 * 3 * % Sigma_row上校% Sigma_col %的背景
函数[f_func, df, num_params] = sym_model_derivatives ()
num_params = 7;
p =符号(“p”[1 num_params]);
信谊f x y A B C tmp_df;
%这里模型定义
一个= cos (p (4)) ^ 2 / (2 * p(6) ^ 2) +罪(p (4)) ^ 2 / (2 * p (5) ^ 2);
B =罪(2 * p (4)) * (1 / p (5) ^ 2 - 1 / p (6) ^ 2) / 4;
C =罪(p (4)) ^ 2 / (2 * p (6) ^ 2) + cos (p (4)) ^ 2 / (2 * p (5) ^ 2);
f = p (1) * exp (- (* (xp (3)) ^ 2 + 2 * B * (xp (3)) * (y p (2)) + C * (y p (2)) ^ 2)) + p (7);
%模型定义的结束
var = [p x y];
为2 = 1:长度(p)
tmp_df = diff (f p (ii));
%转换dfdp matlab函数
df {2} = matlabFunction (tmp_df,“var”一样,var);
结束
% f转换为matlab函数
f_func = matlabFunction (f,“var”一样,var);
结束
% 2 d模型——用于生成数据
函数[f] = model_func (p, x, y, f_func)
p_list = num2cell (p);%,以避免p (1), p (2),…, p (n)参数传递
f = f_func (p_list {:}, x, y);
结束
% Nx2D模型
函数[f] = model_func_ND (c, x, y, f_func)
G_n =大小(x, 2);%数据点的数量
G_k =大小(x, 1);%的数量符合
f = 0 (G_k G_n);
为ff = 1: G_k
:xp = x (ff);%点这适合(n)
:yp = y (ff);
:cp = c (ff);这符合% params(米)
cp_list = num2cell (cp);
% f
:f (ff) = f_func (cp_list {:}, xp, yp);
结束
结束
% Nx2D +动力系统模型
函数(f,动力系统)= model_func_ND_info (c, x, y, f_func)
G_n =大小(x, 2);%数据点的数量
G_k =大小(x, 1);%的数量符合
f = 0 (G_k G_n);
为ff = 1: G_k
:xp = x (ff);%点这适合(n)
:yp = y (ff);
:cp = c (ff);这符合% params(米)
cp_list = num2cell (cp);
% f
:f (ff) = f_func (cp_list {:}, xp, yp);
结束
G_p =元素个数(c) / G_k;%每个适合的参数数量
动力系统= spalloc (G_n * G_k G_p * G_k 0);
结束
% Nx2D +雅可比矩阵
函数f [J] = model_func_ND_J (c, x, y, f_func dfdp_func)
G_n =大小(x, 2);%数据点的数量
G_k =大小(x, 1);%的数量符合
G_p =元素个数(dfdp_func);%每个适合的参数数量
f = 0 (G_k G_n);
J = spalloc (G_n * G_k G_p * G_k G_n * G_k * G_p);
%计算适应值和偏导数为每个健康
为ff = 1: G_k
:xp = x (ff);%点这适合(n)
:yp = y (ff);
:cp = c (ff);这符合% params(米)
cp_list = num2cell (cp);
% f
:f (ff) = f_func (cp_list {:}, xp, yp);
% J
为页= 1:长度(dfdp_func)
df = dfdp_func {pp};%处理函数的偏导数
J (ff: G_k:最终,(pp-1) * G_k + ff) = df (cp_list {:}, xp, yp);
%全面(J)
结束
%大小(J)
结束
结束
% Nx2D +雅可比矩阵相乘
函数(f,动力系统)= model_func_ND_Jm (c xdata ydata、f_func dfdp_func)
全球G_df_func G_P G_n G_k G_xd G_yd;雅可比矩阵乘法函数所需的%
G_n =大小(xdata, 2);%数据点的数量
G_k =大小(xdata, 1);%的数量符合
f = 0 (G_k G_n);
为ff = 1: G_k
:xp = xdata (ff);%点这适合(n)
:yp = ydata (ff);
:cp = c (ff);这符合% params(米)
cp_list = num2cell (cp);
% f
:f (ff) = f_func (cp_list {:}, xp, yp);
结束
%更新全局
G_P = c;%的参数
G_p =元素个数(dfdp_func);%每个适合的参数数量
G_xd = xdata;%指数
G_yd = ydata;
G_df_func = dfdp_func;%偏导数函数处理
%动力系统稀疏矩阵预处理,需要即使雅可比矩阵
%用于繁殖。大小应该等于雅可比矩阵的大小。
动力系统= spalloc (G_n * G_k G_p * G_k G_n * G_k * G_p);%所有零
结束
%计算雅可比矩阵乘法函数。
%注意:雅可比矩阵C是稀疏的。唯一的非零元素
%的衍生品点ξ属于适合j
%对自己的拟合参数。在市场上所处位置
%雅可比矩阵可以通过观察model_func_ND_J雅可比矩阵生成
%使用命令(J)。
函数w = myJ(动力系统,Y,标志)
大小(Y)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
全球G_df_func G_P G_n G_k G_xd G_yd;
P = G_P;
n = G_n;
k = G_k;
xd = G_xd;
码= G_yd;
df_func = G_df_func;
如果国旗> 0
w = Jpositive (Y, df_func, P, n, k, xd,码);
elseif国旗< 0
w = Jnegative (Y, df_func, P, n, k, xd,码);
其他的
w1 = Jpositive (Y, df_func, P, n, k, xd,码);
w = Jnegative (w1、df_func P, n, k, xd,码);
结束
函数= Jpositive (q, df_func P, n, k, xd,码)
%计算C * q(维度NxP Px1 = >资料片)
一个= 0 (n * k, 1);
为ff = 1: k%循环,计算n金额为每个健康
:xp = xd (ff);%点这适合(n)
:yp =码(ff);
:p = p (ff);这符合% params(米)
qp = q (ff:凯西:结束);%的q值符合(m)
油底壳= 0 (1,n);
为页= 1:元素个数(df_func)
df = df_func {pp};% 1 m函数处理
p_list = num2cell (p);
dfp = df (p_list {:}, xp, yp);% n偏导数(可能1值如果导数是一个常量)
油底壳=池+ qp (pp) * dfp;%添加n项为每个参数
结束
(ff:凯西:结束)=池”;% n项
结束
结束
函数= Jnegative (q, df_func P, n, k, xd,码)
%计算C ' * q(维度PxN资料片= > Px1)
= 0 (prod(大小(P)), 1);% Px1
为ff = 1: k%循环,计算3为3适合参数
:xp = xd (ff);%点这适合(n)
:yp =码(ff);
:p = p (ff);这符合% params(米)
qp = q (ff:凯西:结束);%的一部分输入向量问,这符合(n)
为页= 1:元素个数(df_func)
df = df_func {pp};% 1 m函数处理
p_list = num2cell (p);
dfp = df (p_list {:}, xp, yp);% n偏导数
(ff + (pp-1) * k) = (qp的总和。* dfp);
结束
结束
结束
结束
%添加历史如果状态= =“通路”,用于检查收敛
函数停止= myoutput (x, optimvalues状态)
全球历史;
停止= false;
如果isequal(状态,“通路”)
历史=[历史;x];
结束
结束

登录置评。

答案(1)

神骑士
神骑士 2018年11月24日
编辑:神骑士 2018年11月24日
它,而取决于你在做什么。的函数 polyfit 和 interp1 w 工作与 gpuArray 输入。
14日的评论
显示 13年长的评论
布鲁诺陈德良
布鲁诺陈德良 2022年11月10日
编辑:布鲁诺陈德良 2022年11月10日
只是在讨论詹马特的要求在GPU上启用优化功能。
我认为这个功能将mosy有趣的人致力于图像(2 d-3d)处理:去噪,解模糊参数t数量的像素的数量。一些当地运营商的成本函数和非常适合GPU的体系结构。
如果CPU和GPU之间转移数据不需要在每个迭代或每个梯度计算,运行时将会受益。
可能会有一个新的类仍然develipped人员的优化方法,但我认为一个人必须思考更好的利用GPU archiecture TMW优化等等工具。

登录置评。

类别

找到更多的在在MATLAB GPU计算在帮助中心和文件交换

标签

下载188bet金宝搏

社区寻宝

找到宝藏在MATLAB中央,发现社区如何帮助你!

开始狩猎!