Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

Función de multiplicación jacobiana con mínimos cuadrados lineales

Puede resolver un problema de mínimos cuadrados del formulario

minx12Cxd22

tal que lb ≤ x ≤ ub, para los problemas donde es muy grande, tal vez demasiado grande para ser almacenado, mediante el uso de una función de multiplicación jacobiana.C Para esta técnica, utilice el algoritmo.'trust-region-reflective'

Por ejemplo, considere el caso donde es un 2 por matriz basada en una matriz circulante.Cnn Esto significa que las filas de son turnos de un vector de fila.Cv Este ejemplo tiene el vector de fila con elementos del formulario (– 1)vk+1/ :k

= [1, – 1/2, 1/3, – 1/4,..., – 1/],vn

desplazada cíclicamente:

C=[11/21/3...1/n1/n11/2...1/(n1)1/(n1)1/n1...1/(n2)1/21/31/4...111/21/3...1/n1/n11/2...1/(n1)1/(n1)1/n1...1/(n2)1/21/31/4...1].

Este ejemplo de mínimos cuadrados considera el problema donde

= [– 1; 2 ...; – ],dnnn

y las restricciones son – 5 ≤ () ≤ 5 para = 1,...,.xiin

Para lo suficientemente grande, la matriz densa no encaja en la memoria del ordenador. (= 10.000 es demasiado grande en un sistema probado.)nCn

Una función de multiplicación jacobiana tiene la siguiente sintaxis:

w = jmfcn(Jinfo,Y,flag)

es una matriz del mismo tamaño que, usada como preacondicionador.JinfoC Si es demasiado grande para caber en la memoria, debe ser escaso. es un vector o matriz de tamaño por lo que o tiene sentido. indica qué producto se formará:CJinfoYC*YC'*Yflagjmfcn

  • > 0 ⇒flag w = C*Y

  • < 0 ⇒flag w = C'*Y

  • = 0 ⇒flag w = C'*C*Y

Dado que es una matriz tan simple estructurada, es fácil escribir una función de multiplicación jacobiana en términos del vector; es decir, sin formarse.CvC Cada fila de es el producto de una versión desplazada circularmente de veces.C*YvY Se usa para cambiar circularmente.circshiftv

Para calcular, calcular para buscar la primera fila, a continuación, cambiar y calcular la segunda fila, y así sucesivamente.C*Yv*Yv

Para calcular, realice el mismo cálculo, pero utilice una versión desplazada del vector formado a partir de la primera fila de:C'*YtempC'

temp = [fliplr(v),fliplr(v)]; temp = [circshift(temp,1,2),circshift(temp,1,2)]; % Now temp = C'(1,:)

Para calcular, basta con calcular el uso de turnos y, a continuación, calcular los tiempos del resultado utilizando turnos de.C'*C*YC*YvC'fliplr(v)

La función del siguiente código configura el vector y llama al solucionador:dolsqJac3vlsqlin

function [x,resnorm,residual,exitflag,output] = dolsqJac3(n) % r = 1:n-1; % index for making vectors  v(n) = (-1)^(n+1)/n; % allocating the vector v v(r) =( -1).^(r+1)./r;  % Now C should be a 2n-by-n circulant matrix based on v, % but that might be too large to fit into memory.  r = 1:2*n; d(r) = n-r;  Jinfo = [speye(n);speye(n)]; % sparse matrix for preconditioning % This matrix is a required input for the solver; % preconditioning is not really being used in this example  % Pass the vector v so that it does not need to be % computed in the Jacobian multiply function options = optimoptions('lsqlin','Algorithm','trust-region-reflective',...     'JacobianMultiplyFcn',@(Jinfo,Y,flag)lsqcirculant3(Jinfo,Y,flag,v));  lb = -5*ones(1,n); ub = 5*ones(1,n);  [x,resnorm,residual,exitflag,output] = ...     lsqlin(Jinfo,d,[],[],[],[],lb,ub,[],options);

La función de multiplicación jacobiana es la siguiente:lsqcirculant3

function w = lsqcirculant3(Jinfo,Y,flag,v) % This function computes the Jacobian multiply functions % for a 2n-by-n circulant matrix example  if flag > 0     w = Jpositive(Y); elseif flag < 0     w = Jnegative(Y); else     w = Jnegative(Jpositive(Y)); end      function a = Jpositive(q)         % Calculate C*q         temp = v;          a = zeros(size(q)); % allocating the matrix a         a = [a;a]; % the result is twice as tall as the input          for r = 1:size(a,1)             a(r,:) = temp*q; % compute the rth row             temp = circshift(temp,1,2); % shift the circulant         end     end      function a = Jnegative(q)         % Calculate C'*q         temp = fliplr(v);         temp = circshift(temp,1,2); % shift the circulant% the circulant for C'          len = size(q,1)/2; % the returned vector is half as long         % as the input vector         a = zeros(len,size(q,2)); % allocating the matrix a          for r = 1:len             a(r,:) = [temp,temp]*q; % compute the rth row             temp = circshift(temp,1,2); % shift the circulant         end     end end

When = 3000, es una matriz densa de 18 millones elementos.nC Aquí están los resultados de la función para = 3000 en los valores seleccionados de, y la estructura:dolsqJacnxoutput

[x,resnorm,residual,exitflag,output] = dolsqJac3(3000);
 Local minimum possible.  lsqlin stopped because the relative change in function value is less than the function tolerance.
x(1)
ans =     5.0000
x(1500)
ans =    -0.5201
x(3000)
ans =    -5.0000
output
output =     struct with fields:         iterations: 16         algorithm: 'trust-region-reflective'     firstorderopt: 5.9351e-05      cgiterations: 36   constrviolation: []      linearsolver: []           message: 'Local minimum possible.↵↵lsqlin stopped because the relative change in function value is less than the function tolerance.'

Consulte también

|

Temas relacionados