Работа с поворотами Гивенса

Question

Если мы рассмотрим матрицу R размера pxp. Если мы хотим умножить A'RA, где A равно (вращение I + Givens). Здесь I - единичная матрица, а 'обозначает оператор транспонирования.

Мы знаем, что вращение Гивенс является разреженная матрица записывается в виде:

Для выполнения A'RA умножения в MATLAB, мы можем сделать это быстрое внедрение:

%Fast implementation 
    ci = R(:,ik)*(cos(theta))+R(:,jk)*(sin(theta)); % R*A 
    cj = R(:,jk)*(cos(theta)) - R(:,ik)*(sin(theta)); 
    R(:,ik) = ci; 
    R(:,jk) = cj; 

    ri = R(ik,:)*(cos(theta))+R(jk,:)*(sin(theta)); % A'*R*A 
    rj = R(jk,:)*(cos(theta)) - R(ik,:)*(sin(theta)); 
    R(ik,:) = ri; 
    R(jk,:) = rj; 

Но я не понял, как они написали этот код Matlab. Другими словами, я не понимаю, как этот код Matlab применяет умножение A'RA. Пожалуйста, может кто-нибудь помочь мне понять этот код?

Answer 1

Один из возможных источников замешательства заключается в том, что либо знаки в матрице вращения Givens, или сторона, на которой нам нужно транспонировать, в вашем примере неверна. Я предполагаю, что последний: я буду использовать ту же матрицу A, как вы определили, но преобразовать с помощью A*R*A' (изменение транспонирования A эквивалентно принятию угла поворота с противоположным знаком).

Алгоритм относительно прост. Для начала, как комментарии в коде предполагают, преобразование выполняется в два этапа:

Rnew = A * R * A' = A * (R * A') 

Во-первых, мы вычислим R*A'. Для этого представьте матрицу преобразования A = I + M с матрицей вращения Givens M. Формула, которую вы показали в основном, говорит: «Возьмите единичную матрицу, за исключением двух заданных размеров, в которых вы вращаетесь на заданный угол». Вот как полная A матрица выглядит для маленькой задачи (6d матрицы, ik=2, jk=4, как в полной и редкой форме):

Вы можете увидеть, что кроме для (ik, jk) 2d, эта матрица является единичной матрицей, оставляя все остальные измерения неизменными. Таким образом, действие R*A' приведет к R для каждого измерения за исключением для столбцов ik и jk.

В этих двух столбцах результат R*A' является линейной комбинацией R(:,ik) и R(:,jk) с этими тригонометрическими коэффициентами:

[R*A'](:,ik) = R(:,ik)*cos(theta) + R(:,jk)*sin(theta) 
[R*A'](:,jk) = -R(:,ik)*sin(theta) + R(:,jk)*cos(theta) 

, а остальные столбцы остаются неизменными. Если вы посмотрите на код, который вы указали: это именно то, что он делает. Это, по определению, то, что означает R*A', с приведенной выше матрицей A. Все это является следствием того, что матрица A является единичной матрицей, за исключением 2d подпространства.

Следующим шагом является то очень похоже: с помощью этой новой R*A' матрицы умножим из оставил с A.Опять же, эффект по большинству размеров (строк, на этот раз) будет идентичностью, но в рядах ik и jk мы снова получаем линейную комбинацию:

[A*[R*A']](ik,:) = cos(theta)*[R*A'](ik,:) + sin(theta)*[R*A'](jk,:) 
[A*[R*A']](jk,:) = -sin(theta)*[R*A'](ik,:) + cos(theta)*[R*A'](jk,:) 

Отмечая, что код перезаписывает R матрицы с R*A' после первого шага снова становится ясно, что то же самое выполняется в коде «быстрой реализации».

Отказ от ответственности: A' является сопряженным (сопрягаемым транспозисом) в matlab, поэтому вы должны использовать A.' для обозначения транспонирования. Для сложных матриц существует огромная разница, и люди часто забывают использовать правильное транспонирование, когда в конце концов сталкиваются с сложными матрицами.