В этом примере показано, как вычислить 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 строки и любое количество столбцов.
Первый шаг в решении матричного уравнения 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, установите 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-разложение уникально только до признаков строк 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
Используйте типы входных данных с фиксированной точкой для создания эффективного 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>