Алгоритмы, чтобы определить фиксированные точки для действительной матрицы наименьших квадратов решают AX=B

Этот пример показывает алгоритмы что fixed.realQRMatrixSolveFixedpointTypes функционируйте использует, чтобы аналитически определить фиксированные точки для решения действительного матричного уравнения наименьших квадратов 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.realQRMatrixSolveFixedpointTypes функция аналитически вычисляет следующие верхние границы на R, C=QB, и X определить количество целочисленных битов, требуемых избегать переполнения [1,2,3].

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

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

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

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

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

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

Начиная с вычисления svd(A) является более в вычислительном отношении дорогим, чем решение системы уравнений, fixed.realQRMatrixSolveFixedpointTypes функционируйте оценивает нижнюю границу 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+1)Γ(n/2)2m-nΓ(m+12)Γ(m-n+12),m-n+12)

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

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

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

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

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

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

max_abs_A = 1;

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

max_abs_B = 1;

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

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

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

quantizationNoiseStandardDeviation = fixed.realQuantizationNoiseStandardDeviation(precisionBits)
quantizationNoiseStandardDeviation = 1.7206e-08

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

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

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

noiseStandardDeviation = thermalNoiseStandardDeviation;

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

T = fixed.realQRMatrixSolveFixedpointTypes(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: 35
        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 = 17.3205
upperBoundQB = sqrt(m)*max_abs_B
upperBoundQB = 17.3205

Нижняя граница в течение min (svd (A)) для Действительного A

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

estimatedSingularValueLowerBound = fixed.realSingularValueLowerBound(m,n,noiseStandardDeviation)
estimatedSingularValueLowerBound = 0.0371

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

Границы в порядке величины симулированных результатов. Это достаточно, потому что количество битов переводит в логарифмический масштаб относительно области значений значений. Быть в факторе 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 = 17.3205
max(actualMaxR)
ans = 8.3029

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

upperBoundQB
upperBoundQB = 17.3205
max(actualMaxQB)
ans = 2.5707

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

estimatedSingularValueLowerBound
estimatedSingularValueLowerBound = 0.0371
actualSmallestSingularValue = min(singularValues,[],'all')
actualSmallestSingularValue = 0.0420

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

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

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 r e a l 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 r e a l 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.realMatrixSolveUpperBoundX(m,n,max_abs_B,noiseStandardDeviation)
estimated_largest_X = 466.5772
actual_largest_X = max(abs(X_values),[],'all')
actual_largest_X = 44.8056

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

clf
fixed.example.plot.xValueDistribution(m,n,rankA,noiseStandardDeviation,...
    X_values,estimated_largest_X,"real 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*fixed.example.realRandomLowRankMatrix(m,n,rankA);
        % Adding normally distributed random noise makes A non-singular.
        A = A + fixed.example.realNormalRandomArray(0,noiseStandardDeviation,m,n);
        A = quantizenumeric(A,1,A_WordLength,precisionBits);
        B = fixed.example.realUniformRandomArray(-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>

Смотрите также

Функции

Блоки

Похожие темы