fixed.complexQRMatrixSolveFixedpointTypes

Определите фиксированные точки для матричного решения A с комплексным знаком X =B и матричного решения с помощью диагональной загрузки с помощью разложения QR

Описание

пример

T = fixed.realQRMatrixSolveFixedpointTypes(m,n,max_abs_A,max_abs_B,precisionBits) вычисляет фиксированные точки для матричного решения A с комплексным знаком X =B и A λ X =B0, где Aλ=[λIA] и B0=[0B], использование разложение QR. T возвращен как struct с полями, которые задают фиксированные точки для A и B, которые гарантируют, что никакое переполнение не произойдет в QR-алгоритме и X, таким образом, что существует низкая вероятность переполнения.

QR-алгоритм преобразовывает A, оперативный в верхне-треугольный R, и преобразовывает B, оперативный в C =Q'B, где Q R =A является разложением QR A.

пример

T = fixed.realQlessQRMatrixSolveFixedpointTypes(___,noiseStandardDeviation,p_s) задает стандартное отклонение аддитивного случайного шума в A и вероятности, что оценка нижней границы для самого маленького сингулярного значения A больше, чем фактическое самое маленькое сингулярное значение матрицы.

Примеры

свернуть все

Этот пример показывает алгоритмы что fixed.complexQlessQRMatrixSolveFixedpointTypes функционируйте использует, чтобы аналитически определить фиксированные точки для решения комплексного матричного уравнения AAX=B, где A m-n матрица с mn, B n-p, и X n-p.

Обзор

Можно решить матричное уравнение фиксированной точки AAX=B использование разложение QR. Используя последовательность ортогональных преобразований, разложение QR преобразовывает матрицу A оперативный к треугольному верхнему R, где QR=A размер экономики разложение QR. Это уменьшает уравнение до верхне-треугольной системы уравнений RRX=B. Решить для X, вычислить X=R\(R\B) через форварда - и обратная подстановка R в B.

Можно определить соответствующие фиксированные точки для матричного уравнения AAX=B путем выбора дробной длины на основе количества битов точности задан требованиями. fixed.complexQlessQRMatrixSolveFixedpointTypes функция аналитически вычисляет следующие верхние границы на R, и X определить количество целочисленных битов, требуемых избегать переполнения [1,2,3].

Верхняя граница для величины элементов R=QA

max(|R(:)|)mmax(|A(:)|).

Верхняя граница для величины элементов X=(AA)\B

max(|X(:)|)nmax(|B(:)|)min(svd(A))2.

Начиная с вычисления svd(A) является более в вычислительном отношении дорогим, чем решение системы уравнений, fixed.complexQlessQRMatrixSolveFixedpointTypes функционируйте оценивает нижнюю границу min(svd(A)).

Фиксированные точки для решения матричного уравнения (AA)X=B обычно хорошо ограничиваются если количество строк, m, из A очень больше количества столбцов, n т.е. . mn), и A полный ранг. Если A не по сути полный ранг, затем это может быть сделано так путем добавления случайного шума. Случайный шум естественно происходит в физических системах, таких как тепловой шум в радиолокационных системах или системах связи. Если m=n, затем динамический диапазон системы может быть неограниченным, например, в скалярном уравнении x=a2/b и a,b[-1,1]то x может быть произвольно большим если b близко к 0.

Доказательства границ

Свойства и определения векторных и матричных норм

Доказательства границ используют следующие свойства и определения матричных и векторных норм, где Q ортогональная матрица, и v вектор из длины m [6].

||Av||2||A||2||v||2||Q||2=1||v||=max(|v(:)|)||v||||v||2m||v||

Если A m-n матрица и QR=A размер экономики разложение QR A, где Q является ортогональным и m-n и R является верхне-треугольным и n-n, затем сингулярные значения R равны сингулярным значениям A. Если A несингулярно, затем

||R-1||2=||(R)-1||2=1min(svd(R))=1min(svd(A))

Верхняя граница для R = Q'A

Верхняя граница для величины элементов R

max(|R(:)|)mmax(|A(:)|).

Доказательство верхней границы для R = Q'A

jстолбец th R равно R(:,j)=QA(:,j), так

max(|R(:,j)|)=||R(:,j)||||R(:,j)||2=||QA(:,j)||2||Q||2||A(:,j)||2=||A(:,j)||2m||A(:,j)||=mmax(|A(:,j)|)mmax(|A(:)|).

С тех пор max(|R(:,j)|)mmax(|A(:)|) \forall 1jто

max(|R(:)|)mmax(|A(:)|).

Верхняя граница для X = (A'A)\B

Верхняя граница для величины элементов X=(AA)\B

max(|X(:)|)nmax(|B(:)|)min(svd(A))2.

Доказательство верхней границы для X = (A'A)\B

Если A не полный ранг, затем min(svd(A))=0, и если B не равно нулю, затем nmax(|B(:)|)/min(svd(A))2=и таким образом, неравенство верно.

Если AAx=b и QR=A размер экономики разложение QR Aто AAx=RQQRx=RRx=b. Если A полный ранг затем x=R-1((R)-1b)Пусть x=X(:,j) будьте jстолбец th X, и b=B(:,j) будьте jстолбец th Bто

max(|x(:)|)=||x||||x||2=||R-1((R)-1b)||2||R-1||2||(R)-1||2||b||2=(1/min(svd(A))2)||b||2=||b||2/min(svd(A))2n||b||/min(svd(A))2=nmax(|b(:)|)/min(svd(A))2.

С тех пор max(|x(:)|)nmax(|b(:)|)/min(svd(A))2 для всех строк и столбцов B и Xто

max(|X(:)|)nmax(|B(:)|)min(svd(A))2.

Нижняя граница в течение min (svd (A))

Можно оценить нижнюю границу s из min(svd(A))для с комплексным знаком A использование следующей формулы,

s=σN2γ-1(psΓ(m-n+2)2Γ(n)Γ(m+1)Γ(m-n+1)(m-n+1),m-n+1)

где σN стандартное отклонение случайного шума, добавленного к элементам A, 1-ps вероятность это smin(svd(A)), Γ gamma функция, и γ-1обратная неполная гамма функция gammaincinv.

Доказательство найдено в [1]. Это выведено путем интеграции формулы в Лемме 3.4 от [3] и реорганизации терминов.

С тех пор smin(svd(A)) с вероятностью 1-ps, затем вы можете, связал величину элементов X без вычисления svd(A),

max(|X(:)|)nmax(|B(:)|)min(svd(A))2nmax(|B(:)|)s2 с вероятностью 1-ps.

Можно вычислить s использование fixed.complexSingularValueLowerBound функция, которая использует вероятность по умолчанию 5 стандартных отклонений ниже среднего значения, ps=(1+erf(-5/2))/22.866510-7, так вероятность, что предполагаемое направляющееся в самое маленькое сингулярное значение s меньше фактического самого маленького сингулярного значения A 1-ps0.9999997.

Пример

Этот пример запускает симуляцию со многими случайными матрицами и сравнивает аналитические границы с фактическими сингулярными значениями A и фактические самые большие элементы R=QA, и X=(AA)\B.

Задайте системные параметры

Задайте матричные атрибуты и системные параметры для этого примера.

m количество строк в матричном A. В проблеме, такой как beamforming или определение направления, m соответствует количеству отсчетов, которые интегрированы.

m = 300;

n количество столбцов в матричном A и строки в матрицах B и X. В задаче наименьших квадратов, m больше n, и обычно m намного больше, чем n. В проблеме, такой как beamforming или определение направления, n соответствует количеству датчиков.

n = 10;

p количество столбцов в матрицах B и X. Это соответствует одновременному решению системы с p правые стороны.

p = 1;

В этом примере, набор ранг матричного A быть меньше количества столбцов. В проблеме, такой как beamforming или определение направления, rank(A) соответствует количеству сигналов, посягающих на сенсорную матрицу.

rankA = 3;

precisionBits задает количество битов точности, требуемой для матрицы, решают. Установите это значение согласно системным требованиям.

precisionBits = 24;

В этом примере, матрицы с комплексным знаком A и B создаются таким образом, что величина действительных и мнимых частей их элементов меньше чем или равна одному, таким образом, максимальное возможное абсолютное значение любого элемента |1+1i|=2. Ваши собственные системные требования зададут, каковы те значения. Если вы не знаете то, что они, и A и B входные параметры фиксированной точки к системе, затем можно использовать upperbound функция, чтобы определить верхние границы фиксированных точек A и B.

max_abs_A верхняя граница на максимальном элементе массива величины.

max_abs_A = sqrt(2);

max_abs_B верхняя граница на максимальном элементе величины B.

max_abs_B = sqrt(2);

Стандартное отклонение теплового шума является квадратным корнем из степени теплового шума, которая является системным параметром. Хорошо спроектированная система имеет уровень квантования ниже, чем тепловой шум. Здесь, установите thermalNoiseStandardDeviation к эквиваленту -50шумовая мощность дБ.

thermalNoiseStandardDeviation = sqrt(10^(-50/10))
thermalNoiseStandardDeviation = 0.0032

Стандартное отклонение шума от квантования действительных и мнимых частей комплексного сигнала 2-precisionBits/6 [4,5]. Используйте fixed.complexQuantizationNoiseStandardDeviation вычислить это. Смотрите, что это меньше thermalNoiseStandardDeviation.

quantizationNoiseStandardDeviation = fixed.complexQuantizationNoiseStandardDeviation(precisionBits)
quantizationNoiseStandardDeviation = 2.4333e-08

Вычислите фиксированные точки

В этом примере примите что спроектированная системная матрица A не имеет полного ранга (существует меньше сигналов интереса, чем количество столбцов матрицы A), и измеренная системная матрица A имеет аддитивный тепловой шум, который больше, чем шум квантования. Аддитивный шум делает измеренную матрицу A имейте полный ранг.

Набор σnoise=σthermal шум.

noiseStandardDeviation = thermalNoiseStandardDeviation;

Используйте fixed.complexQlessQRMatrixSolveFixedpointTypes вычислить фиксированные точки.

T = fixed.complexQlessQRMatrixSolveFixedpointTypes(m,n,max_abs_A,max_abs_B,...
    precisionBits,noiseStandardDeviation)
T = struct with fields:
    A: [0x0 embedded.fi]
    B: [0x0 embedded.fi]
    X: [0x0 embedded.fi]

T.A тип, вычисленный для преобразования A к R оперативный так, чтобы это не переполнялось.

T.A
ans = 

[]

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 31
        FractionLength: 24

T.B тип, вычисленный для B так, чтобы это не переполнялось.

T.B
ans = 

[]

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 27
        FractionLength: 24

T.X тип, вычисленный для решения X=(AA)\B так, чтобы была низкая вероятность, что это переполняется.

T.X
ans = 

[]

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 40
        FractionLength: 24

Верхняя граница для R

Верхняя граница для R вычисляется с помощью формулы max(|R(:)|)mmax(|A(:)|), где m количество строк матрицы A. Эта верхняя граница используется, чтобы выбрать фиксированную точку с необходимым количеством битов точности, чтобы избежать переполнения в верхней границе.

upperBoundR = sqrt(m)*max_abs_A
upperBoundR = 24.4949

Нижняя граница в течение min (svd (A)) для Комплекса A

Нижняя граница для min(svd(A)) оценивается fixed.complexSingularValueLowerBound функция с помощью вероятности, что оценка s не больше фактического самого маленького сингулярного значения. Вероятность по умолчанию является 5 стандартными отклонениями ниже среднего значения. Можно изменить эту вероятность путем определения его как последнего входного параметра к fixed.complexSingularValueLowerBound функция.

estimatedSingularValueLowerBound = fixed.complexSingularValueLowerBound(m,n,noiseStandardDeviation)
estimatedSingularValueLowerBound = 0.0389

Симулируйте и сравните с вычисленными границами

Границы в порядке величины симулированных результатов. Это достаточно, потому что количество битов переводит в логарифмический масштаб относительно области значений значений. Быть в факторе 10 между 3 и 4 битами. Это - хорошая начальная точка для определения фиксированной точки. Если при запуске симуляцию для большего количества выборок, то более вероятно, что симулированные результаты будут ближе к связанному. Этот пример использует ограниченное количество симуляций, таким образом, не занимает слишком много времени запускаться. Для реальной разработки системы необходимо запустить дополнительные симуляции.

Задайте количество отсчетов, numSamples, по которому можно запустить симуляцию.

numSamples = 1e4;

Запустите симуляцию.

[actualMaxR,singularValues,X_values] = runSimulations(m,n,p,rankA,max_abs_A,max_abs_B,numSamples,...
    noiseStandardDeviation,T);

Вы видите что верхняя граница на R по сравнению с измеренными результатами симуляции максимального значения R по всем запускам в порядке величины.

upperBoundR
upperBoundR = 24.4949
max(actualMaxR)
ans = 9.4990

Наконец, смотрите что предполагаемая нижняя граница min(svd(A)) по сравнению с измеренными результатами симуляции min(svd(A)) по всем запускам также в порядке величины.

estimatedSingularValueLowerBound
estimatedSingularValueLowerBound = 0.0389
actualSmallestSingularValue = min(singularValues,[],'all')
actualSmallestSingularValue = 0.0443

Постройте распределение сингулярных значений по всем запускам симуляции. Распределения самых больших сингулярных значений соответствуют сигналам, которые определяют ранг матрицы. Распределения самых маленьких сингулярных значений соответствуют шуму. Деривация предполагаемого, связанного самого маленького сингулярного значения, использует случайную природу шума.

clf
fixed.example.plot.singularValueDistribution(m,n,rankA,...
    noiseStandardDeviation,singularValues,...
    estimatedSingularValueLowerBound,"complex");

Figure contains an axes object. The axes object with title S i n g u l a r blank v a l u e blank d i s t r i b u t i o n s blank f o r blank 3 0 0 - b y - 1 0 blank c o m p l e x blank m a t r i c e s blank o f blank r a n k blank 3 blank w i t h blank sigma indexOf n o i s e baseline blank = blank 0 . 0 0 3 1 6 contains 20 objects of type line, text.

Увеличьте масштаб к самому маленькому сингулярному значению, чтобы видеть, что связанное предполагаемое близко к нему.

xlim([estimatedSingularValueLowerBound*0.9, max(singularValues(n,:))]);

Figure contains an axes object. The axes object with title S i n g u l a r blank v a l u e blank d i s t r i b u t i o n s blank f o r blank 3 0 0 - b y - 1 0 blank c o m p l e x blank m a t r i c e s blank o f blank r a n k blank 3 blank w i t h blank sigma indexOf n o i s e baseline blank = blank 0 . 0 0 3 1 6 contains 20 objects of type line, text.

Оцените самое большое значение решения, X, и сравните его с самым большим значением X найденный во время запусков симуляции. Оценка в порядке величины фактического значения, которое достаточно для оценки типа данных с фиксированной точкой, потому что это между 3 и 4 битами.

Этот пример использует ограниченное количество запусков симуляции. С запусками дополнительной симуляции фактическое самое большое значение X приблизится к предполагаемому самому большому значению X.

estimated_largest_X = fixed.complexQlessQRMatrixSolveUpperBoundX(m,n,max_abs_B,noiseStandardDeviation)
estimated_largest_X = 9.3348e+03
actual_largest_X = max(abs(X_values),[],'all')
actual_largest_X = 977.7440

Постройте распределение X значений и сравните его с предполагаемой верхней границей для X.

clf
fixed.example.plot.xValueDistribution(m,n,rankA,noiseStandardDeviation,...
    X_values,estimated_largest_X,"complex normally distributed random");

Figure contains an axes object. The axes object with title X blank d i s t r i b u t i o n s blank f o r blank 3 0 0 - b y - 1 0 blank blank m a t r i c e s blank o f blank r a n k blank 3 blank w i t h blank sigma indexOf n o i s e baseline blank = blank 0 . 0 0 3 1 6 contains an object of type line.

Вспомогательные Функции

runSimulations функция создает ряд случайных матриц A и B из данного размера и ранга, квантует их согласно вычисленным типам, вычисляет разложение QR A, и решает уравнение AAX=B. Это возвращает максимальные значения R=QA, сингулярные значения A, и значения X таким образом, их распределения могут быть построены и по сравнению с границами.

function [actualMaxR,singularValues,X_values] = runSimulations(m,n,p,rankA,max_abs_A,max_abs_B,...
        numSamples,noiseStandardDeviation,T)
    precisionBits = T.A.FractionLength;
    A_WordLength = T.A.WordLength;
    B_WordLength = T.B.WordLength;
    actualMaxR = zeros(1,numSamples);
    singularValues = zeros(n,numSamples);
    X_values = zeros(n,numSamples);
    for j = 1:numSamples
        A = (max_abs_A/sqrt(2))*fixed.example.complexRandomLowRankMatrix(m,n,rankA);
        % Adding random noise makes A non-singular.
        A = A + fixed.example.complexNormalRandomArray(0,noiseStandardDeviation,m,n);
        A = quantizenumeric(A,1,A_WordLength,precisionBits);
        B = fixed.example.complexUniformRandomArray(-max_abs_B,max_abs_B,n,p);
        B = quantizenumeric(B,1,B_WordLength,precisionBits);
        [~,R] = qr(A,0);
        X = R\(R'\B);
        actualMaxR(j) = max(abs(R(:)));
        singularValues(:,j) = svd(A);
        X_values(:,j) = X;
    end
end

Ссылки

  1. Томас А. Брайан и Дженна Л. Уоррен. “Системы и Методы для Выбора Расчетного параметра”. Находящийся на рассмотрении патентной заявки. Американская Заявка на патент № 16/947,130. 2020.

  2. Выполните QR-факторизацию Используя CORDIC. Деривация привязанного рост при вычислении QR. MathWorks. 2010. uRL: https://www.mathworks.com/help/fixedpoint/examples/perform-qr-factorization-using-cordic.html.

  3. Цзычжун Чэнь и Джек Дж. Донгарра. “Числа обусловленности Гауссовых Случайных Матриц”. \in: SIAM J. Matrix Anal. Appl. 27.3 (июль 2005), стр 603–620. issn: 0895-4798. doi: 10.1137/040616413. uRL: http://dx.doi.org/10.1137/040616413.

  4. Бернард Видроу. “Исследование Грубого Амплитудного Квантования посредством Теории выборочного метода Найквиста”. \in: Транзакции IRE на Теории 3.4 Схемы (декабрь 1956), стр 266–276.

  5. Бернард Видроу и Иштван Коллар. Шум квантования – ошибка округления в цифровом расчете, обработке сигналов, управлении и коммуникациях. Кембридж, Великобритания: Издательство Кембриджского университета, 2008.

  6. Джин Х. Голуб и Чарльз Ф. ван Лоун. Матричные Расчеты. Второй выпуск. Балтимор: Johns Hopkins University Press, 1989.

Подавите mlint предупреждения в этом файле.

%#ok<*NASGU>
%#ok<*ASGLU> 

Этот пример показывает алгоритмы что fixed.complexQRMatrixSolveFixedpointTypes функционируйте использует, чтобы аналитически определить фиксированные точки для решения комплексного матричного уравнения наименьших квадратов AX=B, где A m-n матрица с mn, B m-p, и X n-p.

Обзор

Можно решить матричное уравнение наименьших квадратов фиксированной точки AX=B использование разложение QR. Используя последовательность ортогональных преобразований, разложение QR преобразовывает матрицу A оперативный к треугольному верхнему R, и преобразовывает матрицу B оперативный к C=QB, где QR=A размер экономики разложение QR. Это уменьшает уравнение до верхне-треугольной системы уравнений RX=C. Решить для X, вычислить X=R\C через заднюю замену R в C.

Можно определить соответствующие фиксированные точки для матричного уравнения наименьших квадратов AX=B путем выбора дробной длины на основе количества битов точности задан требованиями. fixed.complexQRMatrixSolveFixedpointTypes функция аналитически вычисляет следующие верхние границы на R=QA, C=QB, и X определить количество целочисленных битов, требуемых избегать переполнения [1,2,3].

Верхняя граница для величины элементов R=QA

max(|R(:)|)mmax(|A(:)|).

Верхняя граница для величины элементов C=QB

max(|C(:)|)mmax(|B(:)|).

Верхняя граница для величины элементов X=A\B

max(|X(:)|)mmax(|B(:)|)min(svd(A)).

Начиная с вычисления svd(A) является более в вычислительном отношении дорогим, чем решение системы уравнений, fixed.complexQRMatrixSolveFixedpointTypes функционируйте оценивает нижнюю границу min(svd(A)).

Фиксированные точки для решения матричного уравнения AX=B обычно хорошо ограничиваются если количество строк, m, из A очень больше количества столбцов, n т.е. . mn), и A полный ранг. Если A не по сути полный ранг, затем это может быть сделано так путем добавления случайного шума. Случайный шум естественно происходит в физических системах, таких как тепловой шум в радиолокационных системах или системах связи. Если m=n, затем динамический диапазон системы может быть неограниченным, например, в скалярном уравнении x=a/b и a,b[-1,1]то x может быть произвольно большим если b близко к 0.

Доказательства границ

Свойства и определения векторных и матричных норм

Доказательства границ используют следующие свойства и определения матричных и векторных норм, где Q ортогональная матрица, и v вектор из длины m [6].

||Av||2||A||2||v||2||Q||2=1||v||=max(|v(:)|)||v||||v||2m||v||

Если A m-n матрица и QR=A размер экономики разложение QR A, где Q является ортогональным и m-n и R является верхне-треугольным и n-n, затем сингулярные значения R равны сингулярным значениям A. Если A несингулярно, затем

||R-1||2=||(R)-1||2=1min(svd(R))=1min(svd(A))

Верхняя граница для R = Q'A

Верхняя граница для величины элементов R

max(|R(:)|)mmax(|A(:)|).

Доказательство верхней границы для R = Q'A

jстолбец th R равно R(:,j)=QA(:,j), так

max(|R(:,j)|)=||R(:,j)||||R(:,j)||2=||QA(:,j)||2||Q||2||A(:,j)||2=||A(:,j)||2m||A(:,j)||=mmax(|A(:,j)|)mmax(|A(:)|).

С тех пор max(|R(:,j)|)mmax(|A(:)|) \forall 1jто

max(|R(:)|)mmax(|A(:)|).

Верхняя граница для C = Q'B

Верхняя граница для величины элементов C=QB

max(|C(:)|)mmax(|B(:)|).

Доказательство верхней границы для C = Q'B

Доказательство верхней границы для C=QB совпадает с доказательством верхней границы для R=QA путем замены C для R и B для A.

Верхняя граница для X = A\B

Верхняя граница для величины элементов X=A\B

max(|X(:)|)mmax(|B(:)|)min(svd(A)).

Доказательство верхней границы для X = A\B

Если A не полный ранг, затем min(svd(A))=0, и если B не равно нулю, затем mmax(|B(:)|)/min(svd(A))= и таким образом, неравенство верно.

Если A полный ранг, затем x=R-1(Qb)Пусть x=X(:,j) будьте jстолбец th X, и b=B(:,j) будьте jстолбец th Bто

max(|x(:)|)=||x||||x||2=||R-1(Qb)||2||R-1||2||Q||2||b||2=(1/min(svd(A)))1||b||2=||b||2/min(svd(A))m||b||/min(svd(A))=mmax(|b(:)|)/min(svd(A)).

С тех пор max(|x(:)|)mmax(|b(:)|)/min(svd(A)) для всех строк и столбцов B и Xто

max(|X(:)|)mmax(|B(:)|)min(svd(A)).

Нижняя граница в течение min (svd (A))

Можно оценить нижнюю границу s из min(svd(A))для с комплексным знаком A использование следующей формулы,

s=σN2γ-1(psΓ(m-n+2)2Γ(n)Γ(m+1)Γ(m-n+1)(m-n+1),m-n+1)

где σN стандартное отклонение случайного шума, добавленного к элементам A, 1-ps вероятность это smin(svd(A)), Γ gamma функция, и γ-1обратная неполная гамма функция gammaincinv.

Доказательство найдено в [1]. Это выведено путем интеграции формулы в Лемме 3.4 от [3] и реорганизации терминов.

С тех пор smin(svd(A)) с вероятностью 1-ps, затем вы можете, связал величину элементов X без вычисления svd(A),

max(|X(:)|)mmax(|B(:)|)min(svd(A))mmax(|B(:)|)s с вероятностью 1-ps.

Можно вычислить s использование fixed.complexSingularValueLowerBound функция, которая использует вероятность по умолчанию 5 стандартных отклонений ниже среднего значения, ps=(1+erf(-5/2))/22.866510-7, так вероятность, что предполагаемое направляющееся в самое маленькое сингулярное значение s меньше фактического самого маленького сингулярного значения A 1-ps0.9999997.

Пример

Этот пример запускает симуляцию со многими случайными матрицами и сравнивает аналитические границы с фактическими сингулярными значениями A и фактические самые большие элементы R=QA, C=QB, и X=A\B.

Задайте системные параметры

Задайте матричные атрибуты и системные параметры для этого примера.

m количество строк в матрицах A и B. В проблеме, такой как beamforming или определение направления, m соответствует количеству отсчетов, которые интегрированы.

m = 300;

n количество столбцов в матричном A и строки в матричном X. В задаче наименьших квадратов, m больше n, и обычно m намного больше, чем n. В проблеме, такой как beamforming или определение направления, n соответствует количеству датчиков.

n = 10;

p количество столбцов в матрицах B и X. Это соответствует одновременному решению системы с p правые стороны.

p = 1;

В этом примере, набор ранг матричного A быть меньше количества столбцов. В проблеме, такой как beamforming или определение направления, rank(A) соответствует количеству сигналов, посягающих на сенсорную матрицу.

rankA = 3;

precisionBits задает количество битов точности, требуемой для матрицы, решают. Установите это значение согласно системным требованиям.

precisionBits = 24;

В этом примере, матрицы с комплексным знаком A и B создаются таким образом, что величина действительных и мнимых частей их элементов меньше чем или равна одному, таким образом, максимальное возможное абсолютное значение любого элемента |1+1i|=2. Ваши собственные системные требования зададут, каковы те значения. Если вы не знаете то, что они, и A и B входные параметры фиксированной точки к системе, затем можно использовать upperbound функция, чтобы определить верхние границы фиксированных точек A и B.

max_abs_A верхняя граница на максимальном элементе массива величины.

max_abs_A = sqrt(2);

max_abs_B верхняя граница на максимальном элементе величины B.

max_abs_B = sqrt(2);

Стандартное отклонение теплового шума является квадратным корнем из степени теплового шума, которая является системным параметром. Хорошо спроектированная система имеет уровень квантования ниже, чем тепловой шум. Здесь, установите thermalNoiseStandardDeviation к эквиваленту -50шумовая мощность дБ.

thermalNoiseStandardDeviation = sqrt(10^(-50/10))
thermalNoiseStandardDeviation = 0.0032

Стандартное отклонение шума от квантования действительных и мнимых частей комплексного сигнала 2-precisionBits/6 [4,5]. Используйте fixed.complexQuantizationNoiseStandardDeviation функция, чтобы вычислить это. Смотрите, что это меньше thermalNoiseStandardDeviation.

quantizationNoiseStandardDeviation = fixed.complexQuantizationNoiseStandardDeviation(precisionBits)
quantizationNoiseStandardDeviation = 2.4333e-08

Вычислите фиксированные точки

В этом примере примите что спроектированная системная матрица A не имеет полного ранга (существует меньше сигналов интереса, чем количество столбцов матрицы A), и измеренная системная матрица A имеет аддитивный тепловой шум, который больше, чем шум квантования. Аддитивный шум делает измеренную матрицу A имейте полный ранг.

Набор σnoise=σthermal шум.

noiseStandardDeviation = thermalNoiseStandardDeviation;

Используйте fixed.complexQRMatrixSolveFixedpointTypes вычислить фиксированные точки.

T = fixed.complexQRMatrixSolveFixedpointTypes(m,n,max_abs_A,max_abs_B,...
    precisionBits,noiseStandardDeviation)
T = struct with fields:
    A: [0x0 embedded.fi]
    B: [0x0 embedded.fi]
    X: [0x0 embedded.fi]

T.A тип, вычисленный для преобразования A к R оперативный так, чтобы это не переполнялось.

T.A
ans = 

[]

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 31
        FractionLength: 24

T.B тип, вычисленный для преобразования B к QB оперативный так, чтобы это не переполнялось.

T.B
ans = 

[]

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 31
        FractionLength: 24

T.X тип, вычисленный для решения X=A\Bтак, чтобы была низкая вероятность, что это переполняется.

T.X
ans = 

[]

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 36
        FractionLength: 24

Верхние границы для R и C=Q'B

Верхние границы для R и C=QB вычисляются с помощью следующих формул, где m количество строк матриц A и B.

max(|R(:)|)mmax(|A(:)|)

max(|C(:)|)mmax(|B(:)|)

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

upperBoundR = sqrt(m)*max_abs_A
upperBoundR = 24.4949
upperBoundQB = sqrt(m)*max_abs_B
upperBoundQB = 24.4949

Нижняя граница в течение min (svd (A)) для Комплекса A

Нижняя граница для min(svd(A)) оценивается fixed.complexSingularValueLowerBound функция с помощью вероятности, что оценка s не больше фактического самого маленького сингулярного значения. Вероятность по умолчанию является 5 стандартными отклонениями ниже среднего значения. Можно изменить эту вероятность путем определения его как последнего входного параметра к fixed.complexSingularValueLowerBound функция.

estimatedSingularValueLowerBound = fixed.complexSingularValueLowerBound(m,n,noiseStandardDeviation)
estimatedSingularValueLowerBound = 0.0389

Симулируйте и сравните с вычисленными границами

Границы в порядке величины симулированных результатов. Это достаточно, потому что количество битов переводит в логарифмический масштаб относительно области значений значений. Быть в факторе 10 между 3 и 4 битами. Это - хорошая начальная точка для определения фиксированной точки. Если при запуске симуляцию для большего количества выборок, то более вероятно, что симулированные результаты будут ближе к связанному. Этот пример использует ограниченное количество симуляций, таким образом, не занимает слишком много времени запускаться. Для реальной разработки системы необходимо запустить дополнительные симуляции.

Задайте количество отсчетов, numSamples, по которому можно запустить симуляцию.

numSamples = 1e4;

Запустите симуляцию.

[actualMaxR,actualMaxQB,singularValues,X_values] = runSimulations(m,n,p,rankA,max_abs_A,max_abs_B,...
    numSamples,noiseStandardDeviation,T);

Вы видите что верхняя граница на R по сравнению с измеренными результатами симуляции максимального значения R по всем запускам в порядке величины.

upperBoundR
upperBoundR = 24.4949
max(actualMaxR)
ans = 9.6720

Вы видите что верхняя граница на C=QB по сравнению с измеренными результатами симуляции максимального значения C=QB по всем запускам также в порядке величины.

upperBoundQB
upperBoundQB = 24.4949
max(actualMaxQB)
ans = 4.4764

Наконец, смотрите что предполагаемая нижняя граница min(svd(A)) по сравнению с измеренными результатами симуляции min(svd(A)) по всем запускам также в порядке величины.

estimatedSingularValueLowerBound
estimatedSingularValueLowerBound = 0.0389
actualSmallestSingularValue = min(singularValues,[],'all')  
actualSmallestSingularValue = 0.0443

Постройте распределение сингулярных значений по всем запускам симуляции. Распределения самых больших сингулярных значений соответствуют сигналам, которые определяют ранг матрицы. Распределения самых маленьких сингулярных значений соответствуют шуму. Деривация предполагаемого, связанного самого маленького сингулярного значения, использует случайную природу шума.

clf
fixed.example.plot.singularValueDistribution(m,n,rankA,noiseStandardDeviation,...
    singularValues,estimatedSingularValueLowerBound,"complex");

Figure contains an axes object. The axes object with title S i n g u l a r blank v a l u e blank d i s t r i b u t i o n s blank f o r blank 3 0 0 - b y - 1 0 blank c o m p l e x blank m a t r i c e s blank o f blank r a n k blank 3 blank w i t h blank sigma indexOf n o i s e baseline blank = blank 0 . 0 0 3 1 6 contains 20 objects of type line, text.

Увеличьте масштаб к самому маленькому сингулярному значению, чтобы видеть, что связанное предполагаемое близко к нему.

xlim([estimatedSingularValueLowerBound*0.9, max(singularValues(n,:))]);

Figure contains an axes object. The axes object with title S i n g u l a r blank v a l u e blank d i s t r i b u t i o n s blank f o r blank 3 0 0 - b y - 1 0 blank c o m p l e x blank m a t r i c e s blank o f blank r a n k blank 3 blank w i t h blank sigma indexOf n o i s e baseline blank = blank 0 . 0 0 3 1 6 contains 20 objects of type line, text.

Оцените самое большое значение решения, X, и сравните его с самым большим значением X найденный во время запусков симуляции. Оценка в порядке величины фактического значения, которое достаточно для оценки типа данных с фиксированной точкой, потому что это между 3 и 4 битами.

Этот пример использует ограниченное количество запусков симуляции. С запусками дополнительной симуляции фактическое самое большое значение X приблизится к предполагаемому самому большому значению X.

estimated_largest_X = fixed.complexMatrixSolveUpperBoundX(m,n,max_abs_B,noiseStandardDeviation)
estimated_largest_X = 629.3194
actual_largest_X = max(abs(X_values),[],'all')
actual_largest_X = 70.2644

Постройте распределение X значений и сравните его с предполагаемой верхней границей для X.

clf
fixed.example.plot.xValueDistribution(m,n,rankA,noiseStandardDeviation,...
    X_values,estimated_largest_X,"complex normally distributed random");

Figure contains an axes object. The axes object with title X blank d i s t r i b u t i o n s blank f o r blank 3 0 0 - b y - 1 0 blank blank m a t r i c e s blank o f blank r a n k blank 3 blank w i t h blank sigma indexOf n o i s e baseline blank = blank 0 . 0 0 3 1 6 contains an object of type line.

Вспомогательные Функции

runSimulations функция создает ряд случайных матриц A и B из данного размера и ранга, квантует их согласно вычисленным типам, вычисляет разложение QR A, и решает уравнение AX=B. Это возвращает максимальные значения R=QA и C=QB, сингулярные значения A, и значения X таким образом, их распределения могут быть построены и по сравнению с границами.

function [actualMaxR,actualMaxQB,singularValues,X_values] = runSimulations(m,n,p,rankA,max_abs_A,max_abs_B,...
        numSamples,noiseStandardDeviation,T)
    precisionBits = T.A.FractionLength;
    A_WordLength = T.A.WordLength;
    B_WordLength = T.B.WordLength;
    actualMaxR = zeros(1,numSamples);
    actualMaxQB = zeros(1,numSamples);
    singularValues = zeros(n,numSamples);
    X_values = zeros(n,numSamples);
    for j = 1:numSamples    
        A = (max_abs_A/sqrt(2))*fixed.example.complexRandomLowRankMatrix(m,n,rankA);
        % Adding normally distributed random noise makes A non-singular.
        A = A + fixed.example.complexNormalRandomArray(0,noiseStandardDeviation,m,n);
        A = quantizenumeric(A,1,A_WordLength,precisionBits);
        B = fixed.example.complexUniformRandomArray(-max_abs_B,max_abs_B,m,p);
        B = quantizenumeric(B,1,B_WordLength,precisionBits);
        [Q,R] = qr(A,0);
        C = Q'*B;
        X = R\C;
        actualMaxR(j) = max(abs(R(:)));
        actualMaxQB(j) = max(abs(C(:)));
        singularValues(:,j) = svd(A);
        X_values(:,j) = X;
    end
end

Ссылки

  1. Томас А. Брайан и Дженна Л. Уоррен. “Системы и Методы для Выбора Расчетного параметра”. Находящийся на рассмотрении патентной заявки. Американская Заявка на патент № 16/947,130. 2020.

  2. Выполните QR-факторизацию Используя CORDIC. Деривация привязанного рост при вычислении QR. MathWorks. 2010. uRL: https://www.mathworks.com/help/fixedpoint/examples/perform-qr-factorization-using-cordic.html.

  3. Цзычжун Чэнь и Джек Дж. Донгарра. “Числа обусловленности Гауссовых Случайных Матриц”. \in: SIAM J. Matrix Anal. Appl. 27.3 (июль 2005), стр 603–620. issn: 0895-4798. doi: 10.1137/040616413. uRL: http://dx.doi.org/10.1137/040616413.

  4. Бернард Видроу. “Исследование Грубого Амплитудного Квантования посредством Теории выборочного метода Найквиста”. \in: Транзакции IRE на Теории 3.4 Схемы (декабрь 1956), стр 266–276.

  5. Бернард Видроу и Иштван Коллар. Шум квантования – ошибка округления в цифровом расчете, обработке сигналов, управлении и коммуникациях. Кембридж, Великобритания: Издательство Кембриджского университета, 2008.

  6. Джин Х. Голуб и Чарльз Ф. ван Лоун. Матричные Расчеты. Второй выпуск. Балтимор: Johns Hopkins University Press, 1989.

Подавите mlint предупреждения в этом файле.

%#ok<*NASGU>
%#ok<*ASGLU>

В этом примере показано, как использовать fixed.complexQRMatrixSolveFixedpointTypes функционируйте, чтобы аналитически определить фиксированные точки для решения комплексного матричного уравнения наименьших квадратов AX=B, где A m-n матрица с mn, B m-p, и X n-p.

Фиксированные точки для решения матричного уравнения AX=B хорошо ограничены если количество строк, m, из A очень больше количества столбцов, n т.е. . mn), и A полный ранг. Если A не по сути полный ранг, затем это может быть сделано так путем добавления случайного шума. Случайный шум естественно происходит в физических системах, таких как тепловой шум в радиолокационных системах или системах связи. Если m=n, затем динамический диапазон системы может быть неограниченным, например, в скалярном уравнении x=a/b и a,b[-1,1]то x может быть произвольно большим если b близко к 0.

Задайте системные параметры

Задайте матричные атрибуты и системные параметры для этого примера.

m количество строк в матрицах A и B. В проблеме, такой как beamforming или определение направления, m соответствует количеству отсчетов, которые интегрированы.

m = 300;

n количество столбцов в матричном A и строки в матричном X. В задаче наименьших квадратов, m больше n, и обычно m намного больше, чем n. В проблеме, такой как beamforming или определение направления, n соответствует количеству датчиков.

n = 10;

p количество столбцов в матрицах B и X. Это соответствует одновременному решению системы с p правые стороны.

p = 1;

В этом примере, набор ранг матричного A быть меньше количества столбцов. В проблеме, такой как beamforming или определение направления, rank(A) соответствует количеству сигналов, посягающих на сенсорную матрицу.

rankA = 3;

precisionBits задает количество битов точности, требуемой для матрицы, решают. Установите это значение согласно системным требованиям.

precisionBits = 24;

В этом примере, матрицы с комплексным знаком A и B создаются таким образом, что величина действительных и мнимых частей их элементов меньше чем или равна одному, таким образом, максимальное возможное абсолютное значение любого элемента |1+1i|=2. Ваши собственные системные требования зададут, каковы те значения. Если вы не знаете то, что они, и A и B входные параметры фиксированной точки к системе, затем можно использовать upperbound функция, чтобы определить верхние границы фиксированных точек A и B.

max_abs_A верхняя граница на максимальном элементе массива величины.

max_abs_A = sqrt(2);  

max_abs_B верхняя граница на максимальном элементе величины B.

max_abs_B = sqrt(2);

Стандартное отклонение теплового шума является квадратным корнем из степени теплового шума, которая является системным параметром. Хорошо спроектированная система имеет уровень квантования ниже, чем тепловой шум. Здесь, установите thermalNoiseStandardDeviation к эквиваленту -50шумовая мощность дБ.

thermalNoiseStandardDeviation = sqrt(10^(-50/10))
thermalNoiseStandardDeviation = 0.0032

Стандартное отклонение шума квантования является функцией необходимого количества битов точности. Используйте fixed.complexQuantizationNoiseStandardDeviation вычислить это. Смотрите, что это меньше thermalNoiseStandardDeviation.

quantizationNoiseStandardDeviation = fixed.complexQuantizationNoiseStandardDeviation(precisionBits)
quantizationNoiseStandardDeviation = 2.4333e-08

Вычислите фиксированные точки

В этом примере примите что спроектированная системная матрица A не имеет полного ранга (существует меньше сигналов интереса, чем количество столбцов матрицы A), и измеренная системная матрица A имеет аддитивный тепловой шум, который больше, чем шум квантования. Аддитивный шум делает измеренную матрицу A имейте полный ранг.

Набор σnoise=σthermal шум.

noiseStandardDeviation = thermalNoiseStandardDeviation;

Используйте fixed.complexQRMatrixSolveFixedpointTypes вычислить фиксированные точки.

T = fixed.complexQRMatrixSolveFixedpointTypes(m,n,max_abs_A,max_abs_B,...
    precisionBits,noiseStandardDeviation)
T = struct with fields:
    A: [0x0 embedded.fi]
    B: [0x0 embedded.fi]
    X: [0x0 embedded.fi]

T.A тип, вычисленный для преобразования A к R=QA оперативный так, чтобы это не переполнялось.

T.A
ans = 

[]

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 31
        FractionLength: 24

T.B тип, вычисленный для преобразования B к C=QB оперативный так, чтобы это не переполнялось.

T.B
ans = 

[]

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 31
        FractionLength: 24

T.X тип, вычисленный для решения X=A\Bтак, чтобы была низкая вероятность, что это переполняется.

T.X
ans = 

[]

          DataTypeMode: Fixed-point: binary point scaling
            Signedness: Signed
            WordLength: 36
        FractionLength: 24

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

Создайте случайные матрицы A и B таким образом, что B находится в области значений A, и rankA=rank(A). Добавьте случайный шум измерения в A который заставит его стать полным рангом, но это будет также влиять на решение так, чтобы B только близко к области значений A.

rng('default');
[A,B] = fixed.example.complexRandomLeastSquaresMatrices(m,n,p,rankA);
A = A + fixed.example.complexNormalRandomArray(0,noiseStandardDeviation,m,n);

Бросьте входные параметры к типам, определенным fixed.complexQRMatrixSolveFixedpointTypes. Квантование к фиксированной точке эквивалентно добавлению случайного шума.

A = cast(A,'like',T.A);
B = cast(B,'like',T.B);

Ускорьте fixed.qrMatrixSolve функция при помощи fiaccel сгенерировать исполняемый файл MATLAB (MEX) функция.

fiaccel fixed.qrMatrixSolve -args {A,B,T.X} -o qrComplexMatrixSolve_mex

Задайте выходной тип T.X и вычислите фиксированную точку X=A\B использование метода QR.

X = qrComplexMatrixSolve_mex(A,B,T.X);

Вычислите относительную погрешность, чтобы проверить точность выхода.

relative_error = norm(double(A*X - B))/norm(double(B))
relative_error = 0.0056

Подавите mlint предупреждения в этом файле.

%#ok<*NASGU>
%#ok<*ASGLU>

Входные параметры

свернуть все

Количество строк в A и B в виде положительного скаляра с целочисленным знаком.

Типы данных: double

Количество столбцов в A в виде положительного скаляра с целочисленным знаком.

Типы данных: double

Максимум абсолютного значения A в виде скаляра.

Пример: max(abs(A(:)))

Типы данных: double

Максимум абсолютного значения B в виде скаляра.

Пример: max(abs(B(:)))

Типы данных: double

Необходимое количество битов точности ввода и вывода в виде положительного скаляра с целочисленным знаком.

Типы данных: double

Стандартное отклонение аддитивного случайного шума в A в виде скаляра.

Если noiseStandardDeviation не задан, затем значением по умолчанию является стандартное отклонение шума квантования с комплексным знаком σq=(2precisionBits)/(6), который вычисляется fixed.complexQuantizationNoiseStandardDeviation.

Типы данных: double

Вероятность, что оценка нижней границы s больше, чем фактическое самое маленькое сингулярное значение матрицы в виде скаляра. Использование fixed.complexSingularValueLowerBound оценить самое маленькое сингулярное значение, s, A. Если p_s не задан, значение по умолчанию ps=(1/2)(1+erf(5/2))3107 который является 5 стандартными отклонениями ниже среднего значения, таким образом, вероятностью, что предполагаемое направляющееся в самое маленькое сингулярное значение меньше фактического самого маленького сингулярного значения, является 1-ps ≈ 0.9999997.

Типы данных: double

Выходные аргументы

свернуть все

Фиксированные точки для A, B, и X, возвратились как struct. Struct T имеет поля T.A, T.B, и T.X. Эти поля содержат fi объекты, которые задают фиксированные точки для

  • A и B, которые не гарантируют переполнения, произойдут в QR-алгоритме.

    QR-алгоритм преобразовывает A, оперативный в верхне-треугольный R, и преобразовывает B, оперативный в C =Q'B, где Q R =A является разложением QR A.

  • X, таким образом, что существует низкая вероятность переполнения.

Советы

Использование fixed.realQRMatrixSolveFixedpointTypes вычислить фиксированные точки для входных параметров этих функций и блоков.

Алгоритмы

T.A и T.B вычисляются с помощью fixed.qrFixedpointTypes. Количество целочисленных битов, требуемых предотвратить переполнение, выведено из следующих границ на росте R и C =Q'B [1]. Необходимое количество целочисленных битов добавляется к количеству битов точности, precisionBits, из входа, плюс один для знакового бита, плюс один бит для промежуточного усиления CORDIC приблизительно 1,6468 [2].

Элементы R ограничены в величине

max(|R(:)|)mmax(|A(:)|).

Элементы C =Q'B ограничены в величине

max(|C(:)|)mmax(|B(:)|).

T.X вычисляется путем ограничения выхода, X, в решении методом наименьших квадратов A X =B использование следующей формулы [3] [4].

Элементы X =R \(Q 'B) ограничены в величине

max(|X(:)|)mmax(|B(:)|)min(svd(A)).

Вычисление сингулярного разложения, чтобы вывести вышеупомянутое привязало X, является более в вычислительном отношении дорогим, чем целая матрица решает, таким образом, fixed.complexSingularValueLowerBound функция используется, чтобы оценить привязанный min(svd(A)).

Ссылки

[2] Voler, Джек Э. "Тригонометрический вычислительный метод CORDIC". Транзакции IRE на электронно-вычислительных машинах EC-8 (1959): 330-334.

[3] Брайан, Томас А. и Дженна Л. Уоррен. "Системы и методы для выбора расчетного параметра". Американская заявка на патент № 16/947, 130. 2020.

[4] Чен, Цзычжун и Джек Дж. Донгарра. "Числа обусловленности Гауссовых Случайных Матриц". SIAM Journal на Анализе матрицы и Приложения 27, № 3 (июль 2005): 603-620.

Введенный в R2021b
Для просмотра документации необходимо авторизоваться на сайте