Обзор

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

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

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

$$\max (|R(:)|)\le \sqrt{m}\max (|A(:)|)$$.

Верхняя граница для величины элементов $$C=Q^{\prime }B$$

$$\max (|C(:)|)\le \sqrt{m}\max (|B(:)|)$$.

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

$$\max (|X(:)|)\le \frac{\sqrt{m}\max (|B(:)|)}{\min (\text{svd}(A))}$$.

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

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

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

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

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

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

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

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

$$\max (|R(:)|)\le \sqrt{m}\max (|A(:)|)$$.

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

$$j$$столбец th $$R$$ равно $$R(:,j)=Q^{\prime }A(:,j)$$, так

С тех пор $$\max (|R(:,j)|)\le \sqrt{m}\max (|A(:)|)$$ \forall $$1\le j$$то

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

Верхняя граница для величины элементов $$C=Q^{\prime }B$$

$$\max (|C(:)|)\le \sqrt{m}\max (|B(:)|)$$.

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

Доказательство верхней границы для $$C=Q^{\prime }B$$ совпадает с доказательством верхней границы для $$R=Q^{\prime }A$$ путем замены $$C$$ для $$R$$ и $$B$$ для $$A$$.

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

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

$$\max (|X(:)|)\le \frac{\sqrt{m}\max (|B(:)|)}{\min (\text{svd}(A))}$$.

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

Если $$A$$ не полный ранг, затем $$\min (\text{svd}(A))=0$$, и если $$B$$ не равно нулю, затем $$\sqrt{m}\max (|B(:)|)/\min (\text{svd}(A))=\infty$$ и таким образом, неравенство верно.

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

С тех пор $$\max (|x(:)|)\le \sqrt{m}\max (|b(:)|)/\min (\text{svd}(A))$$ для всех строк и столбцов $$B$$ и $$X$$то

$$\max (|X(:)|)\le \frac{\sqrt{m}\max (|B(:)|)}{\min (\text{svd}(A))}$$.

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

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

где $${\sigma }_N$$ стандартное отклонение случайного шума, добавленного к элементам $$A$$, $$1-p_s$$ вероятность это $$s\le \min (\text{svd}(A))$$, $$\Gamma$$ gamma функция, и $${\gamma }^{-1}$$обратная неполная гамма функция gammaincinv.

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

С тех пор $$s\le \min (\text{svd}(A))$$ с вероятностью $$1-p_s$$, затем вы можете, связал величину элементов $$X$$ без вычисления $$\text{svd}(A)$$,

$$\max (|X(:)|)\le \frac{\sqrt{m}\max (|B(:)|)}{\min (\text{svd}(A))}\le \frac{\sqrt{m}\max (|B(:)|)}{s}$$ с вероятностью $$1-p_s$$.

Можно вычислить $$s$$ использование fixed.realSingularValueLowerBound функция, которая использует вероятность по умолчанию 5 стандартных отклонений ниже среднего значения $$p_s=(1+\text{erf}(-5/\sqrt{2}))/2\approx 2.8665\cdot 10^{-7}$$, так вероятность, что предполагаемое направляющееся в самое маленькое сингулярное значение $$s$$ меньше фактического самого маленького сингулярного значения $$A$$ $$1-p_s\approx 0.9999997$$.

Пример

Этот пример запускает симуляцию со многими случайными матрицами и сравнивает аналитические границы с фактическими сингулярными значениями $$A$$ и фактические самые большие элементы $$R=Q^{\prime }A$$, $$C=Q^{\prime }B$$, и $$X=A\backslash 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 или определение направления, $$\text{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^{-\text{precisionBits}}/\sqrt{12}$$ [4,5]. Используйте fixed.realQuantizationNoiseStandardDeviation функция, чтобы вычислить это. Смотрите, что это меньше thermalNoiseStandardDeviation.

quantizationNoiseStandardDeviation = fixed.realQuantizationNoiseStandardDeviation(precisionBits)

quantizationNoiseStandardDeviation = 1.7206e-08

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

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

Набор $${\sigma }_{\text{noise}}={\sigma }_{\text{thermal}\text{шум}}$$.

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 тип, вычисленный для преобразования $\mathit{A}$ к $\mathit{R}$ оперативный так, чтобы это не переполнялось.

T.A

ans = 

[]

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

T.B тип, вычисленный для преобразования $\mathit{B}$ к ${\mathit{Q}}^{\prime }\mathit{B}$ оперативный так, чтобы это не переполнялось.

T.B

ans = 

[]

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

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

T.X

ans = 

[]

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

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

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

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

upperBoundR = sqrt(m)*max_abs_A

upperBoundR = 17.3205

upperBoundQB = sqrt(m)*max_abs_B

upperBoundQB = 17.3205

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

Нижняя граница для $$\min (\text{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=Q^{\prime }B$$ по сравнению с измеренными результатами симуляции максимального значения $$C=Q^{\prime }B$$ по всем запускам также в порядке величины.

upperBoundQB

upperBoundQB = 17.3205

max(actualMaxQB)

ans = 2.5707

Наконец, смотрите что предполагаемая нижняя граница $$\min (\text{svd}(A))$$ по сравнению с измеренными результатами симуляции $$\min (\text{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");

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

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

Оцените самое большое значение решения, 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");

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

runSimulations функция создает ряд случайных матриц $$A$$ и $$B$$ из данного размера и ранга, квантует их согласно вычисленным типам, вычисляет разложение QR $$A$$, и решает уравнение $$AX=B$$. Это возвращает максимальные значения $$R=Q^{\prime }A$$ и $$C=Q^{\prime }B$$, сингулярные значения $$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

Ссылки

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

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

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

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

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 или определение направления, $$\text{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

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

quantizationNoiseStandardDeviation = fixed.realQuantizationNoiseStandardDeviation(precisionBits)

quantizationNoiseStandardDeviation = 1.7206e-08

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

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

Набор $${\sigma }_{\text{noise}}={\sigma }_{\text{thermal}\text{шум}}$$.

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 тип, вычисленный для преобразования $\mathit{A}$ к $\mathit{R}={\mathit{Q}}^{\prime }\mathit{A}$ оперативный так, чтобы это не переполнялось.

T.A

ans = 

[]

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

T.B тип, вычисленный для преобразования $\mathit{B}$ к ${\mathit{C}=\mathit{Q}}^{\prime }\mathit{B}$ оперативный так, чтобы это не переполнялось.

T.B

ans = 

[]

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

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

T.X

ans = 

[]

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

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

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

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

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

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 qrRealMatrixSolve_mex

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

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

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

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

relative_error = 0.0063

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

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