Вычислить взвешенное суммирование матричной мощности (матричный полином) в Matlab

Имея nxn-матрицу A_k и nx1-вектор x, есть ли какой-нибудь умный способ вычислить

введите здесь описание изображения

с помощью Матлаба? x_i — элементы вектора x, поэтому J — сумма матриц. До сих пор я использовал цикл for, но мне было интересно, есть ли более разумный способ.


matlab matrix-multiplication polynomials
person cholo14    schedule 04.05.2017    source источник
comment
Что означает A_k^i? А (я, к)? Можете ли вы привести небольшой пример Matlab?   -  person m7913d    schedule 04.05.2017
comment
Извините, A_k — это просто имя матрицы, поэтому A_k^i — степень A_k для i.   -  person cholo14    schedule 04.05.2017
comment
Так что же такое J, это матрица или число? Являются ли {A_k} формой последовательности матриц?   -  person Miriam Farber    schedule 04.05.2017
comment
И x_i - это вектор nx1 или i-й элемент x?   -  person m7913d    schedule 04.05.2017

Ответы (1)


arrow_upward
1
arrow_downward

Короткий ответ: вы можете использовать встроенную функцию Matlab polyvalm для вычисления матричного полинома следующим образом:

x = x(end:-1:1); % flip the order of the elements
x(end+1) = 0; % append 0
J = polyvalm(x, A);

Длинный ответ: внутри Matlab используется цикл. Таким образом, вы не так много выиграли или работали еще хуже, если оптимизировали собственную реализацию (см. мою функцию calcJ_loopOptimised):

% construct random input
n = 100;
A = rand(n);
x = rand(n, 1);

% calculate the result using different methods
Jbuiltin = calcJ_builtin(A, x);
Jloop = calcJ_loop(A, x);
JloopOptimised = calcJ_loopOptimised(A, x);

% check if the functions are mathematically equivalent (should be in the order of `eps`)
relativeError1 = max(max(abs(Jbuiltin - Jloop)))/max(max(Jbuiltin))
relativeError2 = max(max(abs(Jloop - JloopOptimised)))/max(max(Jloop))

% measure the execution time
t_loopOptimised = timeit(@() calcJ_loopOptimised(A, x))
t_builtin = timeit(@() calcJ_builtin(A, x))
t_loop = timeit(@() calcJ_loop(A, x))

% check if builtin function is faster
builtinFaster = t_builtin < t_loopOptimised

% calculate J using Matlab builtin function
function J = calcJ_builtin(A, x)
  x = x(end:-1:1);
  x(end+1) = 0;
  J = polyvalm(x, A);
end

% naive loop implementation
function J = calcJ_loop(A, x)
  n = size(A, 1);
  J = zeros(n,n);
  for i=1:n
    J = J + A^i * x(i);
  end
end

% optimised loop implementation (cache result of matrix power)
function J = calcJ_loopOptimised(A, x)
  n = size(A, 1);
  J = zeros(n,n);
  A_ = eye(n);
  for i=1:n
    A_ = A_*A;
    J = J + A_ * x(i);
  end
end

Для n=100 я получаю следующее:

t_loopOptimised = 0.0077
t_builtin       = 0.0084
t_loop          = 0.0295

Для n=5 я получаю следующее:

t_loopOptimised = 7.4425e-06
t_builtin       = 4.7399e-05
t_loop          = 1.0496e-04

Обратите внимание, что мои тайминги несколько колеблются между разными запусками, но оптимизированный цикл почти всегда быстрее (до 6 раз для маленьких n), чем встроенная функция.

person m7913d    schedule 04.05.2017