QR алгоритм нахождения собственных значений

@cmath · Регистрация: 11.08.2012

Студворк — интернет-сервис помощи студентам

Доброго времени суток. Такая у меня маленькая проблемка:
Реализую алгоритм, указанный в заголовке. Как разложить матрицу на ортогональную*верхнетреугольную дошел, не дошел, когда выходить из цикла
A1=R0*Q0
A2=R1*Q1 (A1=Q1*R1)
...
Вообще, как я понял, надо выйти из цикла, как только матрица будет An=треуголного вида (с некоторой точностью).
Но: если у матрицы в спектре есть комплекснозначные собственные числа, то треугольного вида не будет

. Матрица будет квазитреугольная, т.е. будут блоки на главной диагонали. Смотрел на Учитесь.ру инфу по этому вопросу, но не сильно помогло. Подскажите плиз, каким критерием пользоваться для выхода из цикла в общем случае.

@011 · 10.12.2015, 02:19

cmath, скажите пожалуйста, доделали вы реализацию и не могли бы поделиться?

@cmath · 10.12.2015, 02:26 **[ТС]**

Сообщение от 011

cmath, скажите пожалуйста, доделали вы реализацию и не могли бы поделиться?

Нет, не доделал. Я тогда только условия себе на бумаге выписал - как закончить, но до конкретной реализации дело не дошло. С преподом договорились, что остановимся на матрицах, в спектре которых нет комплексных чисел. Впрочем сложность состоит конкретно в проверке "квазитреугольности", Если матрица имеет четную размерность, то тут все просто, а если нечетную, есть варианты размещения блоков. У меня был алгоритм для общего случая, но затерялся в архивах. Если найду и не забуду - выложу. Реализацию с "треугольным" случаем выложу часиков через 10, когда с работы приду.

@011 · 10.12.2015, 08:04

cmath, а у меня матрица нечетной размерности (7).

Скажите, а как вы вначале алгоритма раскладываете А? (Т.е. как вы находите Q, R? A = QR) Преобразованиями Хаусхолдера?
Еще интересует такой вопрос: все-таки когда нужно прекратить итерации? Видел такой случай: продолжают, пока поддиагональные элементы становятся меньшими заданной точности.

Не по теме:

Спасибо, что откликнулись :)

@cmath · 10.12.2015, 08:16 **[ТС]**

Сообщение от 011

Преобразованиями Хаусхолдера?

Они самые, родимые.

Сообщение от 011

все-таки когда нужно прекратить итерации? Видел такой случай: продолжают, пока поддиагональные элементы становятся меньшими заданной точности.

Это я и реализовал. Но заковыка: если в спектре есть числа с мнимой частью, то данный критерий не сработает и программа зацикливается. Поэтому я её тестировал только на матрицах, про которые железно известно, что все собственные числа - действительные - например, на симметрических матрицах (думаю объяснять этот момент не нужно).

@011 · 10.12.2015, 10:51

cmath, я бы очень хотел посмотреть на Вашу реализацию

@cmath · 10.12.2015, 12:10 **[ТС]**

Готово. Смотрим архив там две проги (проекты VC++2008). В одной - нахождение собственных чисел матрицы 3х3, в другой просто разложение QR. Делал давненько, когда еще только начинал C++ учить, все работает, но ничему там не удивляйтесь

@011 · 10.12.2015, 14:23

cmath, Здорово, благодарю! Только мне надо будет перенести на размерность 5. И в связи с этим есть вопрос: в проекте "Алгоритм QR" как нужно изменить подсчет a, b, участвующих в определении точности?

C++
1
2
3
4
        a = pow(A[1][0]*A[1][0]+A[2][0]*A[2][0], 0.5);
        b = fabs(A[2][1]);
        if(a < eps && b < eps) ex = false;
        else ex = true;

@cmath · 10.12.2015, 15:28 **[ТС]**

Там производится суммирование квадратов поддиагональных элементов по столбцу и извлекается корень. Отсюда вывод...

Добавлено через 4 минуты
Да, кстати, для случая с размерностью =5 там будет 4 числа. Можно, канешь, обобщить алгоритм на матрицы произвольной размерности, но у меня сейчас нет времени этим заниматься.

@011 · 10.12.2015, 20:24

cmath,
А почему А[2][1] не участвовал в сумме квадратов под корнем?

Правильно ли я понимаю, что в данном случае (все СЗ -- действительные) можно суммировать все поддиагональные элементы? (вернее находить корень из сумм их квадратов).
Это повысит точность оценки процесса?

И еще вопрос: блоки возникают, когда у нас есть комплексные СЗ? Выше вы говорили об алгоритме на проверку квазитреугольности матрицы: удалось ли вам его найти?
У моей матрицы все СЗ -- действительные, но просто интересно, как этот алгоритм выглядит

@cmath · 11.12.2015, 06:22 **[ТС]**

Сообщение от 011

Это повысит точность оценки процесса?

Да, можно вообще сравнивать корень из суммы квадратов всех поддиагональных элементов, в таком случае при выполнении неравенства (корень меньше эпсилон) все эти числа должны быть еще меньше, чем если сравнивать с эпсилон максимальный по модулю элемент. Я брал корни сумм квадратов элементов по столбцам, вследствие этого A[2][1] в первую сумму и не попал, зато попал во вторую (число b), написал модуль, т.к. корень из квадрата числа равен его модулю (поскольку других слагаемых нет).

Добавлено через 1 час 4 минуты

Сообщение от 011

алгоритме на проверку квазитреугольности матрицы

Нужно искать блоки на главной диагонали. Есть у матрицы главная диагональ, если рассматривать следующую (что под ней) диагональ, то уже под ней должны элементы быть меньше эпсилона на самой этой диагонали для каждого элемента, большего, чем эпсилон, два "соседа" (тот, что выше на строку и тот, что ниже) должны быть меньше эпсилон. Таким способом можно "подцепить" блок "за уголок" и вытянуть все остальные его элементы, причем найденный угловой элемент блока будет левым нижним углом, что нужно учесть при получении индексов других элементов блока. Удачи.

Добавлено через 6 минут
P.S. Последний алгоритм можно применить к матрице с размерностью любой четности.

@011 · 11.12.2015, 12:51

cmath, не могли бы Вы помочь с переносом вашего кода на размерность 5х5?

Я попробовал добавить "итерации разложения" по аналогии с первой и второй. Но что-то не так...

Измененный (5х5)

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
#include <windows.h>
#include <cmath>
#include <iostream>
using namespace std;
 
int sign(double arg)
{
    if(arg > 0) return 1;
    if(arg == 0) return 0;
    if(arg < 0) return -1;
}
 
const int n = 5;
 
int main()
{
    setlocale(LC_ALL, "Russian");
 
    freopen("in.txt", "r", stdin);
 
    double A[n][n] = {0}; // Матрица, для которой будет определён спектр
    double B[n][n] = {0}; // Вспомогательная матрица
    double H[n][n] = {0}; // Матрица Хаусхолдера
    double Q[n][n] = {0}; // Матрица Q разложения
    double V[n]    = {0}; // Вектор V
    double Num     =   0; // Скалярное произведение вектора самого на себя
    // Заполнение матриц
    for(int i = 0; i < n; i++)
    {
        Q[i][i] = 1;
        H[i][i] = 1;
    }
    //cout << "Введите матрицу A 5x5:" << endl;
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < n; j++)
        {
            //cout << "a[" << i+1 << "][" << j+1 << "] = "; 
            cin >> A[i][j]; 
        }
    }
    cout << "A 5x5:" << endl;
    for(int i = 0; i < n; i++)
    {
        for(int j = 0; j < n; j++)
        {
            cout << A[i][j];
            if(j != 2) cout << " ";
            else cout << endl;
        }
    }
    double eps;
    double a, b;
    //cout << "Введите точность:"; 
    cin >> eps;
    bool ex = true;
    // Вход в цикл
    while(ex)
    {
        // первая итерация разложения
        // нахождение вектора V
        V[0] = A[0][0] + sign(A[0][0])*pow(A[0][0] * A[0][0] + A[1][0] * A[1][0] + A[2][0] * A[2][0] + A[3][0] * A[3][0] + A[4][0] * A[4][0], 0.5);
        V[1] = A[1][0];
        V[2] = A[2][0];
        V[3] = A[3][0];
        V[4] = A[4][0];
        // квадрат евклидовой нормы V
        Num = V[0]*V[0] + V[1]*V[1] + V[2]*V[2] + V[3]*V[3] + V[4]*V[4];
        // получение матрицы Хаусхолдера
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                H[i][j] -= 2*V[i]*V[j]/Num;
                Q[i][j] = H[i][j];
            }
        }
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                for(int k = 0; k < n; k++)
                {
                    B[i][j] += H[i][k]*A[k][j];
                }
            }
        }
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                A[i][j] = B[i][j];
                B[i][j] = 0;
                if(i != j) H[i][j] = 0;
                else H[i][j] = 1;
            }
        }
        // вторая итерация разложения
        V[0] = 0; 
        V[1] = A[1][1] + sign(A[1][1])*pow(A[1][1] * A[1][1] + A[2][1] * A[2][1] + A[3][1]*A[3][1] + A[4][1]*A[4][1], 0.5);
        V[2] = A[2][1];
        V[3] = A[3][1];
        V[4] = A[4][1];
        Num = V[0] * V[0] + V[1] * V[1] + V[2] * V[2] + V[3] * V[3] + V[4] * V[4];
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                H[i][j] -= 2*V[i]*V[j]/Num;
                B[i][j] = Q[i][j]; // сохраняем в B первоначальное значение Q
                Q[i][j] = 0;
            }
        }
        // третья итерация разложения
        V[0] = 0;
        V[1] = 0;
        V[2] = A[2][2] + sign(A[2][2])*pow(A[2][2] * A[2][2] + A[3][2] * A[3][2] + A[4][2] * A[4][2], 0.5);
        V[3] = A[3][1];
        V[4] = A[4][1];
        Num = V[0] * V[0] + V[1] * V[1] + V[2] * V[2] + V[3] * V[3] + V[4] * V[4];
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                H[i][j] -= 2 * V[i] * V[j] / Num;
                B[i][j] = Q[i][j]; // сохраняем в B первоначальное значение Q
                Q[i][j] = 0;
            }
        }
        // четвертая итерация разложения
        V[0] = 0;
        V[1] = 0;
        V[2] = 0;
        V[3] = A[3][3] + sign(A[3][3])*pow(A[3][3] * A[3][3] + A[3][4] * A[3][4], 0.5);
        V[4] = A[4][1];
        Num = V[0] * V[0] + V[1] * V[1] + V[2] * V[2] + V[3] * V[3] + V[4] * V[4];
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                H[i][j] -= 2 * V[i] * V[j] / Num;
                B[i][j] = Q[i][j]; // сохраняем в B первоначальное значение Q
                Q[i][j] = 0;
            }
        }
 
        // получение Q (завершающий этап)
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                for(int k = 0; k < n; k++)
                {
                    Q[i][j] += B[i][k]*H[k][j];
                }
            }
        }
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                B[i][j] = 0;
            }
        }
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                for(int k = 0; k < n; k++)
                {
                    B[i][j] += H[i][k]*A[k][j];
                }
            }
        }
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                A[i][j] = B[i][j];
                B[i][j] = 0;
            }
        }
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                for(int k = 0; k < n; k++)
                {
                    B[i][j] += A[i][k]*Q[k][j];
                }
            }
        }
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                A[i][j] = B[i][j];
                B[i][j] = 0;
            }
        }
        for(int i = 0; i < n; i++)
        {
            for(int j = 0; j < n; j++)
            {
                if(i != j) H[i][j] = 0;
                else H[i][j] = 1;
            }
        }
        a = pow(A[1][0] * A[1][0] + A[2][0] * A[2][0] + A[2][1] * A[2][1] + A[3][0] * A[3][0] + A[3][1] * A[3][1] + A[3][2] * A[3][2] + A[4][0] * A[4][0] + A[4][1] * A[4][1] + A[4][2] * A[4][2] + A[4][3] * A[4][3], 0.5);
        //b = fabs(A[2][1]);
        if(a < eps/* && b < eps*/) ex = true;
        else ex = false;
    }
    cout << "Собственные значения:";
    //cout << " L[1] = " << A[0][0] << "; L[2] = " << A[1][1] << "; L[3] = " << A[2][2] << endl;
    for (int i = 0; i < n; ++i)
        cout << A[i][i] << endl;
    //cout << "a = " << a << endl; cout << "b = " << b << endl;
    system("pause");
    return 0;
}

in.txt

Да, я понимаю, что double не хранит столько знаков. Но, если получится на double, то на длинную арифметику я перенесу.

Code
1
2
3
4
5
6
7
94.7823411408704000 29.7947350625651000 -27.0238605223580000    50.8020105795783000 1.2295806527721800
29.7947350625651000 44.5735589546075000 -39.3874753663386000    21.2721127377941000 -0.2745837863773420
-27.0238605223580000    -39.3874753663386000    41.5292744513661000 -36.7017978223677000    -9.3405335929985800
50.8020105795783000 21.2721127377941000 -36.7017978223677000    81.5067295079263000 26.2417822118307000
1.2295806527721800  -0.2745837863773420 -9.3405335929985800 26.2417822118307000 26.3240078779145000
 
0.5

Исходный (3х3)

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
#include "main.h"
 
using namespace std;
 
int sign(double arg)
{
    if(arg > 0) return 1;
    if(arg == 0) return 0;
    if(arg < 0) return -1;
}
 
const int n = 5;
 
int main()
{
    setlocale(LC_ALL, "Russian");
    double A[3][3] = {0}; // Матрица, для которой будет определён спектр
    double B[3][3] = {0}; // Вспомогательная матрица
    double H[3][3] = {0}; // Матрица Хаусхолдера
    double Q[3][3] = {0}; // Матрица Q разложения
    double V[3]    = {0}; // Вектор V
    double Num     =   0; // Скалярное произведение вектора самого на себя
    // Заполнение матриц
    for(int i = 0; i < 3; i++)
    {
        Q[i][i] = 1;
        H[i][i] = 1;
    }
    cout << "Введите матрицу A 3x3:" << endl;
    for(int i = 0; i < 3; i++)
    {
        for(int j = 0; j < 3; j++)
        {
            cout << "a[" << i+1 << "][" << j+1 << "] = "; cin >> A[i][j]; 
        }
    }
    cout << "A 3x3:" << endl;
    for(int i = 0; i < 3; i++)
    {
        for(int j = 0; j < 3; j++)
        {
            cout << A[i][j];
            if(j != 2) cout << " ";
            else cout << endl;
        }
    }
    double eps;
    double a, b;
    cout << "Введите точность:"; cin >> eps;
    bool ex = true;
    // Вход в цикл
    while(ex)
    {
        // первая итерация разложения
        // нахождение вектора V
        V[0] = A[0][0]+sign(A[0][0])*pow(A[0][0]*A[0][0]+A[1][0]*A[1][0]+A[2][0]*A[2][0],0.5);
        V[1] = A[1][0];
        V[2] = A[2][0];
        // квадрат евклидовой нормы V
        Num = V[0]*V[0] + V[1]*V[1] + V[2]*V[2];
        // получение матрицы Хаусхолдера
        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                H[i][j] -= 2*V[i]*V[j]/Num;
                Q[i][j] = H[i][j];
            }
        }
        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                for(int k = 0; k < 3; k++)
                {
                    B[i][j] += H[i][k]*A[k][j];
                }
            }
        }
        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                A[i][j] = B[i][j];
                B[i][j] = 0;
                if(i != j) H[i][j] = 0;
                else H[i][j] = 1;
            }
        }
        // вторая итерация разложения
        V[1] = A[1][1]+sign(A[1][1])*pow(A[1][1]*A[1][1]+A[2][1]*A[2][1], 0.5);
        V[0] = 0;
        V[2] = A[2][1];
        Num = V[1]*V[1] + V[2]*V[2];
        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                H[i][j] -= 2*V[i]*V[j]/Num;
                B[i][j] = Q[i][j]; // сохраняем в B первоначальное значение Q
                Q[i][j] = 0;
            }
        }
        // получение Q (завершающий этап)
        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                for(int k = 0; k < 3; k++)
                {
                    Q[i][j] += B[i][k]*H[k][j];
                }
            }
        }
        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                B[i][j] = 0;
            }
        }
        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                for(int k = 0; k < 3; k++)
                {
                    B[i][j] += H[i][k]*A[k][j];
                }
            }
        }
        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                A[i][j] = B[i][j];
                B[i][j] = 0;
            }
        }
        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                for(int k = 0; k < 3; k++)
                {
                    B[i][j] += A[i][k]*Q[k][j];
                }
            }
        }
        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                A[i][j] = B[i][j];
                B[i][j] = 0;
            }
        }
        for(int i = 0; i < 3; i++)
        {
            for(int j = 0; j < 3; j++)
            {
                if(i != j) H[i][j] = 0;
                else H[i][j] = 1;
            }
        }
        a = pow(A[1][0]*A[1][0]+A[2][0]*A[2][0], 0.5);
        b = fabs(A[2][1]);
        if(a < eps && b < eps) ex = false;
        else ex = true;
    }
    cout << "Собственные значения:";
    cout << " L[1] = " << A[0][0] << "; L[2] = " << A[1][1] << "; L[3] = " << A[2][2] << endl;
    cout << "a = " << a << endl; cout << "b = " << b << endl;
    system("pause");
    return 0;
}

@cmath · 11.12.2015, 13:43 **[ТС]**

Не по теме:

Лень-матушка...

Второй проект открывали? Хотя бы ради интереса. Там, меж прочим, разложение QR (без поиска СЗ) для общего случая (тоже не идеал, но все-таки), можно было бы использовать тот код и не писать //четвертая итерация и т.д. Допуская массу ошибок. Сам подход, который я тогда использовал вследствие отсутствия опыта и навыков - топорный, грабля короче, на которую вы сейчас наступаете. Я дал вам код второго проекта не просто так - проанализируйте его, там есть все, что нужно.
P.S. Уже обладая некоторым опытом и навыками я бы на вашем месте сначала бы написал класс для квадратных матриц и для векторов, переопределил бы для них операции сложения, умножения, код бы понятнее стал на порядок.

iifat · 11.12.2015, 16:54

Сообщение от cmath

Подскажите плиз, каким критерием пользоваться для выхода из цикла в общем случае

Как я понял, матрица должна сходиться к двухдиагональной. На главной собственные числа, диагональю выше — единицы с нулями. Нули в местах перехода от одного собственного числа к другому.

Новые блоги и статьи Все статьи Все блоги /
Кому нужен AOT? DevAlt 26.03.2026 Решил сделать простой ланчер Написал заготовку: dotnet new console --aot -o UrlHandler var items = args. Split(":"); var tag = items; var id = items; var executable = args;. . .	Отправка уведомления на почту при изменении наименования справочника Maks 24.03.2026 Программная отправка письма электронной почты на примере изменения наименования типового справочника "Склады" в конфигурации БП3. Перед реализацией необходимо выполнить настройку системной учетной. . .	модель ЗдравоСохранения 5. Меньше увольнений- больше дохода! anaschu 24.03.2026 Теперь система здравосохранения уменьшает количество увольнений. 9TO2GP2bpX4 a42b81fb172ffc12ca589c7898261ccb/ https:/ / rutube. ru/ video/ a42b81fb172ffc12ca589c7898261ccb/ Слева синяя линия -. . .	Midnight Chicago Blues kumehtar 24.03.2026 Такой Midnight Chicago Blues, знаешь?. . Когда вечерние улицы становятся ночными, а ты не можешь уснуть. Ты идёшь в любимый старый бар, и бармен наливает тебе виски. Ты смотришь на пролетающие. . .
SDL3 для Desktop (MinGW): Вывод текста со шрифтом TTF с помощью библиотеки SDL3_ttf на Си и C++ 8Observer8 24.03.2026 Содержание блога Финальные проекты на Си и на C++: finish-text-sdl3-c. zip finish-text-sdl3-cpp. zip	Жизнь в неопределённости kumehtar 23.03.2026 Жизнь — это постоянное существование в неопределённости. Например, даже если у тебя есть список дел, невозможно дойти до точки, где всё окончательно завершено и больше ничего не осталось. В принципе,. . .	Модель здравоСохранения: работники работают быстрее после её введения. anaschu 23.03.2026 geJalZw1fLo Корпорация до введения программа здравоохранения имела много невыполненных работниками заданий, после введения программы количество заданий выросло. Но на выплатах по больничным это. . .	Контроль уникальности заводского номера Maks 23.03.2026 Алгоритм контроля уникальности заводского (или серийного) номера на примере нетипового документа выдачи шин для спецтехники с табличной частью, разработанного в конфигурации КА2. Номеклатура. . .

@cmath 2525 / 1751 / 152 Регистрация: 11.08.2012 Сообщений: 3,349

	QR алгоритм нахождения собственных значений 23.10.2012, 16:11. Показов 11251. Ответов 13 Метки нет (Все метки) Доброго времени суток. Такая у меня маленькая проблемка: Реализую алгоритм, указанный в заголовке. Как разложить матрицу на ортогональнуюверхнетреугольную дошел, не дошел, когда выходить из цикла A1=R0Q0 A2=R1Q1 (A1=Q1R1) ... Вообще, как я понял, надо выйти из цикла, как только матрица будет An=треуголного вида (с некоторой точностью). Но: если у матрицы в спектре есть комплекснозначные собственные числа, то треугольного вида не будет. Матрица будет квазитреугольная, т.е. будут блоки на главной диагонали. Смотрел на Учитесь.ру инфу по этому вопросу, но не сильно помогло. Подскажите плиз, каким критерием пользоваться для выхода из цикла в общем случае. 0

@011 10 / 10 / 1 Регистрация: 28.11.2013 Сообщений: 153
	10.12.2015, 02:19
	cmath, скажите пожалуйста, доделали вы реализацию и не могли бы поделиться? 0

@011 10 / 10 / 1 Регистрация: 28.11.2013 Сообщений: 153
	10.12.2015, 08:04
	cmath, а у меня матрица нечетной размерности (7). Скажите, а как вы вначале алгоритма раскладываете А? (Т.е. как вы находите Q, R? A = QR) Преобразованиями Хаусхолдера? Еще интересует такой вопрос: все-таки когда нужно прекратить итерации? Видел такой случай: продолжают, пока поддиагональные элементы становятся меньшими заданной точности. Не по теме: Спасибо, что откликнулись :) 0

@011 10 / 10 / 1 Регистрация: 28.11.2013 Сообщений: 153
	10.12.2015, 10:51
	cmath, я бы очень хотел посмотреть на Вашу реализацию 0

@cmath 2525 / 1751 / 152 Регистрация: 11.08.2012 Сообщений: 3,349
	10.12.2015, 15:28 [ТС]
	Там производится суммирование квадратов поддиагональных элементов по столбцу и извлекается корень. Отсюда вывод... Добавлено через 4 минуты Да, кстати, для случая с размерностью =5 там будет 4 числа. Можно, канешь, обобщить алгоритм на матрицы произвольной размерности, но у меня сейчас нет времени этим заниматься. 1

@011 10 / 10 / 1 Регистрация: 28.11.2013 Сообщений: 153
	10.12.2015, 20:24
	cmath, А почему А[2][1] не участвовал в сумме квадратов под корнем? Правильно ли я понимаю, что в данном случае (все СЗ -- действительные) можно суммировать все поддиагональные элементы? (вернее находить корень из сумм их квадратов). Это повысит точность оценки процесса? И еще вопрос: блоки возникают, когда у нас есть комплексные СЗ? Выше вы говорили об алгоритме на проверку квазитреугольности матрицы: удалось ли вам его найти? У моей матрицы все СЗ -- действительные, но просто интересно, как этот алгоритм выглядит 0

@cmath 2525 / 1751 / 152 Регистрация: 11.08.2012 Сообщений: 3,349
	11.12.2015, 13:43 [ТС]
	Не по теме: Лень-матушка... Второй проект открывали? Хотя бы ради интереса. Там, меж прочим, разложение QR (без поиска СЗ) для общего случая (тоже не идеал, но все-таки), можно было бы использовать тот код и не писать //четвертая итерация и т.д. Допуская массу ошибок. Сам подход, который я тогда использовал вследствие отсутствия опыта и навыков - топорный, грабля короче, на которую вы сейчас наступаете. Я дал вам код второго проекта не просто так - проанализируйте его, там есть все, что нужно. P.S. Уже обладая некоторым опытом и навыками я бы на вашем месте сначала бы написал класс для квадратных матриц и для векторов, переопределил бы для них операции сложения, умножения, код бы понятнее стал на порядок. 0