Реализуйте аппаратно-эффективное QR-разложение с помощью CORDIC в Systolic Array

В этом примере показано, как вычислить QR-разложение матриц с помощью аппаратно эффективного кода MATLAB ® в Simulink ®.

Чтобы решить систему уравнений или вычислить решение методом наименьших квадратов матричного уравнения AX = B с помощью QR-разложения, вычислите R и Q 'B, где QR = A, и RX = Q' B. R является верхней треугольной матрицей, а Q - ортогональной матрицей. Если вы просто хотите Q и R, тогда установите B, чтобы быть единичной матрицей.

В этом примере R вычисляется из матрицы А путем применения преобразований Гивенса с помощью алгоритма CORDIC. C = Q 'B вычисляется из матрицы B путем применения тех же преобразований Гивенса. Алгоритм использует только итерационные сдвиги и сложения, чтобы выполнить эти расчеты.

Для получения дополнительной информации об алгоритме, используемом в этом примере, смотрите Выполните QR-факторизацию Используя CORDIC

Обзор

Модель Simulink, используемая в этом примере, является:

fxpdemo_real_4x4_systolic_array_QR_model

Чтобы ввести свои собственные матрицы входа A и B, откройте параметры блоков соответствующих включенных блоков подсистемы налево. После симуляции модель возвращает вычисленные выходные матрицы C = Q 'B и R в рабочую область. Можно задать количество итераций CORDIC в блочных параметрах подсистемы транспонирования (Q) * B, R 4x4 Real CORDIC Systolic Array. Если входы являются фиксированной точкой, то количество итераций CORDIC должно быть меньше размера слова. Точность расчета улучшает один бит для каждой итерации, вплоть до размера слова входов.

Эта модель будет работать с типами данных с фиксированной точкой, двойной точкой и одним типом.

Чтобы увидеть, как алгоритм выполняет факторизацию, смотрите под маской подсистемы транспонирования (Q) * B, R 4x4 Real CORDIC Systolic Array. Аннотации указывают, какие строки матрицы работают для обнуления поддиагональных элементов. Систематический массив настраивается для матрицы А 4 на 4, но может быть расширен до любого размера путем следования тому же шаблону. Эта реализация работает только с действительными входными матрицами.

Чтобы увидеть код MATLAB в блоке MATLAB Function, который выполняет преобразования Givens с помощью CORDIC, продолжите искать под масками блоков.

В этом примере количество строк и столбцов A должно быть 4. Матрица B должна иметь 4 строки и любое количество столбцов.

Используйте QR, чтобы решить матричное уравнение Ax = B

Первый шаг в решении матричного уравнения AX = B - вычислить RX = Q 'B, где R является верхнетреугольной, Q является ортогональным и Q * R = A.

Следующие входы являются типами с плавающей точностью с двойной точностью, поэтому установите количество итераций равным 52, что является количеством бит в мантиссе double.

format
NumberOfCORDICIterations = 52;
A = 2*rand(4,4)-1;
B = 2*rand(4,4)-1;

Симулируйте модель, чтобы вычислить R и C = Q 'B.

sim fxpdemo_real_4x4_systolic_array_QR_model
R
C
R =

    1.5149   -0.0519    1.7292   -0.3224
         0    0.9593   -0.0259   -0.0879
         0         0    0.2565    1.0888
         0         0         0   -0.6429


C =

    0.5942   -0.2382    0.0676   -0.9370
   -0.8887    0.6146   -0.5758    0.3051
    0.1725    0.7339    0.5409    0.5374
    0.8540    1.1078   -0.2183   -0.5620

Проверьте, что обратная замена R и C = Q 'B дает те же результаты, что и MATLAB.

X = R\C
X_should_be = A\B
X =

   -7.1245  -12.1131   -0.6637    1.4236
   -0.8779    0.7572   -0.5511    0.3545
    6.3113   10.1759    0.6673   -1.6155
   -1.3283   -1.7231    0.3396    0.8741


X_should_be =

   -7.1245  -12.1131   -0.6637    1.4236
   -0.8779    0.7572   -0.5511    0.3545
    6.3113   10.1759    0.6673   -1.6155
   -1.3283   -1.7231    0.3396    0.8741

Норма различия между встроенным MATLAB и Решением CORDIC должна быть небольшой.

norm(X - X_should_be)
ans =

   3.6171e-14

Вычисление Q и R

Чтобы вычислить Q и R, установите B равной матрице тождеств. Когда B равна единичной матрице, то Q = C '.

NumberOfCORDICIterations = 52;
A = 2*rand(4,4)-1;
B = eye(size(A,1),'like',A);
sim fxpdemo_real_4x4_systolic_array_QR_model

Q = C';

Теоретическое QR-разложение является QR = = A, поэтому различие между вычисленными QR и A должна быть небольшой.

norm(Q*R - A)
ans =

   2.2861e-15

QR не является уникальным

QR-разложение уникально только до признаков строк R и столбцов Q. Вы можете сделать уникальное QR-разложение, сделав диагональные элементы R всеми положительными.

D = diag(sign(diag(R)));
Qunique = Q*D
Runique = D*R
Qunique =

   -0.3086    0.1224   -0.1033   -0.9376
   -0.6277   -0.7636   -0.0952    0.1174
   -0.5573    0.3930    0.7146    0.1559
    0.4474   -0.4975    0.6852   -0.2877


Runique =

    1.4459   -0.8090    0.1547    0.3977
         0    1.1441    0.0809   -0.2494
         0         0    0.8193    0.1894
         0         0         0    0.4836

Затем можно сравнить вычисленный QR из модели с функцией builtin MATLAB qr.

[Q0,R0] = qr(A);
D0 = diag(sign(diag(R0)));
Q0 = Q0*D0
R0 = D0*R0
Q0 =

   -0.3086    0.1224   -0.1033   -0.9376
   -0.6277   -0.7636   -0.0952    0.1174
   -0.5573    0.3930    0.7146    0.1559
    0.4474   -0.4975    0.6852   -0.2877


R0 =

    1.4459   -0.8090    0.1547    0.3977
         0    1.1441    0.0809   -0.2494
         0         0    0.8193    0.1894
         0         0         0    0.4836

Используйте Fixed-Point для аппаратно-эффективной реализации

Используйте типы входных данных с фиксированной точкой для создания эффективного HDL-кода для устройств ASIC и FPGA.

Дополнительные сведения о том, как выбрать типы данных с фиксированной точкой, которые не будут переполнены, см. в примере Выполнение QR-факторизации с использованием CORDIC.

Можно запустить много тестовых воздействий через модель, сделав A и B 3-мерные массивы.

n_test_inputs=100;

В следующем разделе заданы случайные входы для матриц А и B, которые масштабируются между -1 и 1, поэтому установите размер слова с фиксированной точкой в 18 битах и длину дроби в 14 битах, чтобы позволить рост в QR-факторизации и промежуточных расчетах в алгоритме CORDIC.

word_length = 18;
fraction_length = 14;

Наиболее точным числом итераций CORDIC является размер слова минус один. Если для количества итераций CORDIC задано значение меньше word_length - 1, то задержка и тактовые такты для следующего сигнала готовности будут короче, но это будет менее точно. Количество итераций CORDIC не должно быть больше, поскольку сгенерированный код не поддерживает сдвиги, большие, чем размер слова типа с фиксированной точкой.

NumberOfCORDICIterations = word_length - 1
NumberOfCORDICIterations =

    17

Случайные тестовые воздействия объединены так, что в момент времени k, входы являются A (:,:, k) и B (:,:, k). Каждый элемент массива и B является равномерной случайной переменной между -1 и + 1.

A = 2*rand(4,4,n_test_inputs)-1;

Выберите B, чтобы быть единичной матрицей, поэтому Q = C '.

B = eye(4);
B = repmat(B,1,1,n_test_inputs);

Приведите A к фиксированной точке и приведите B как A.

A = fi(A,1,word_length,fraction_length);
B = cast(B,'like',A);

Симулируйте модель

sim fxpdemo_real_4x4_systolic_array_QR_model

Вычислите и постройте график ошибок

norm_error = zeros(1,size(R,3));
for k = 1:size(R,3)
    Q_times_R_minus_A = double(C(:,:,k))'*double(R(:,:,k)) - double(A(:,:,k));
    norm_error(k) = norm(Q_times_R_minus_A);
end

Ошибки должны быть порядка 10 ^ -3.

clf
plot(norm_error,'o-')
grid on
title('norm(QR - A)')

%#ok<*NASGU,*NOPTS>
Для просмотра документации необходимо авторизоваться на сайте