A block-recursive approach to unitary matrix decomposition
Singular value decomposition is one of the fundamental tools of numerical linear algebra and is widely used in data processing, optimization, computer vision, and machine learning. The classical approach to its computation is based on the Householder algorithm, which reduces an arbitrary matrix to a...
Saved in:
| Date: | 2026 |
|---|---|
| Main Authors: | , |
| Format: | Article |
| Language: | Ukrainian |
| Published: |
PROBLEMS IN PROGRAMMING
2026
|
| Subjects: | |
| Online Access: | https://pp.isofts.kiev.ua/index.php/ojs1/article/view/895 |
| Tags: |
Add Tag
No Tags, Be the first to tag this record!
|
| Journal Title: | Problems in programming |
| Download file: |  |
Institution
Problems in programming| _version_ | 1863311599691890688 |
|---|---|
| author | Malaschonok, G.I. Sukharskyi, S.S. |
| author_facet | Malaschonok, G.I. Sukharskyi, S.S. |
| author_sort | Malaschonok, G.I. |
| baseUrl_str | https://pp.isofts.kiev.ua/index.php/ojs1/oai |
| collection | OJS |
| datestamp_date | 2026-04-23T22:26:13Z |
| description | Singular value decomposition is one of the fundamental tools of numerical linear algebra and is widely used in data processing, optimization, computer vision, and machine learning. The classical approach to its computation is based on the Householder algorithm, which reduces an arbitrary matrix to a tridiagonal or bidiagonal form followed by the computation of singular values. However, as the matrix size increases, significant difficulties arise in efficiently parallelizing such algorithms, which limits their performance on modern high-performance computing systems. One promising approach to overcoming these limitations is the use of block-recursive algorithms. Such algo rithms allow the original problem to be decomposed into independent block subproblems that can be processed in parallel, providing a high degree of scalability on multiprocessor systems and distributed-memory clusters. This paper proposes a new block-recursive algorithm for the orthogonal (or unitary) decomposition of symmetric matrices based on QR and QP decompositions. The algorithm constructs orthogonal (unitary) factors and pro duces a band matrix, reducing the computational complexity of subsequent stages of spectral decomposition. A complexity analysis is presented showing that the asymptotic order of the algorithm is determined by the complexity of matrix multiplication.Experimental performance studies were carried out in a parallel computing environment using GPU acceleration for matrix operations. The obtained results demonstrate the efficiency of the proposed algorithm as the problem size increases and confirm the agreement between the experimental com plexity and the theoretical estimates. The reconstruction error of the matrix remains at the level of machine precision for double-precision floating-point numbers, which indicates the numerical stability of the proposed approach.Problems in programming 2026; 1: 102-111 |
| first_indexed | 2026-04-24T01:00:16Z |
| format | Article |
| fulltext |
102
Аналітика даних
© Г.І. Малашонок, С.С. Сухарський, 2026
ISSN 1727-4907. Проблеми програмування. 2026. №1
УДК 519.61:004.272 https://doi.org/10.15407/pp2026.01.102
Г. І. Малашонок, С. С. Сухарський
БЛОКОВО-РЕКУРСИВНИЙ ПІДХІД ДО УНІТАРНОГО
РОЗКЛАДУ МАТРИЦЬ
Сингулярне розкладання матриці є одним з основних інструментів чисельної лінійної алгебри та широко
застосовується у задачах обробки даних, оптимізації, комп’ютерного зору та машинного навчання. Кла-
сичний підхід до його обчислення базується на алгоритмі Хаусхолдера, який дозволяє звести довільну
матрицю до тридіагональної або бідіагональної форми з подальшим обчисленням сингулярних значень.
Однак зі зростанням розмірності матриць виникають суттєві труднощі з ефективним розпаралелюванням
таких алгоритмів, що обмежує їхню продуктивність на сучасних високопродуктивних обчислювальних
системах.
Одним із перспективних підходів до подолання цих обмежень є використання блоково-рекурсивних ал-
горитмів. Такі алгоритми дозволяють розділяти початкову задачу на незалежні блокові підзадачі, що мо-
жуть виконуватися паралельно, забезпечуючи високий ступінь масштабованості на багатопроцесорних
системах та кластерах із розподіленою пам’яттю.
У даній роботі запропоновано новий алгоритм ортогонального (унітарного) блоково-рекурсивного роз-
кладу симетричних матриць, який базується на використанні QR- та QP-розкладів. Алгоритм дозволяє
побудувати ортогональні (унітарні) множники та отримати стрічкову матрицю, що зменшує обчислюва-
льну складність подальших етапів спектрального розкладу.
Отримано оцінку обчислювальної складності алгоритму та показано, що її асимптотичний порядок ви-
значається складністю множення матриць. Експериментальні дослідження продуктивності проведено у
середовищі паралельних обчислень із використанням прискорення матричних операцій на графічному
процесорі. Отримані результати демонструють ефективність запропонованого алгоритму при зростанні
розміру задачі та підтверджують узгодженість експериментальної складності з теоретичними оцінками.
Значення похибки реконструкції матриці знаходяться на рівні машинної точності для чисел подвійної
точності, що свідчить про чисельну стабільність запропонованого підходу.
Ключові слова: унітарний розклад матриці, блоково-рекурсивний алгоритм, паралельні обчислення, гра-
фічний процесор
G. I. Malaschonok, S. S. Sukharskyi
A BLOCK-RECURSIVE APPROACH TO UNITARY
MATRIX DECOMPOSITION
Singular value decomposition is one of the fundamental tools of numerical linear algebra and is widely used in
data processing, optimization, computer vision, and machine learning. The classical approach to its computation
is based on the Householder algorithm, which reduces an arbitrary matrix to a tridiagonal or bidiagonal form
followed by the computation of singular values. However, as the matrix size increases, significant difficulties
arise in efficiently parallelizing such algorithms, which limits their performance on modern high-performance
computing systems.
One promising approach to overcoming these limitations is the use of block-recursive algorithms. Such algo-
rithms allow the original problem to be decomposed into independent block subproblems that can be processed
in parallel, providing a high degree of scalability on multiprocessor systems and distributed-memory clusters.
This paper proposes a new block-recursive algorithm for the orthogonal (or unitary) decomposition of symmetric
matrices based on QR and QP decompositions. The algorithm constructs orthogonal (unitary) factors and pro-
duces a band matrix, reducing the computational complexity of subsequent stages of spectral decomposition.
A complexity analysis is presented showing that the asymptotic order of the algorithm is determined by the
complexity of matrix multiplication.Experimental performance studies were carried out in a parallel computing
environment using GPU acceleration for matrix operations. The obtained results demonstrate the efficiency of
the proposed algorithm as the problem size increases and confirm the agreement between the experimental com-
plexity and the theoretical estimates. The reconstruction error of the matrix remains at the level of machine
precision for double-precision floating-point numbers, which indicates the numerical stability of the proposed
approach.
Keywords: unitary matrix decomposition, block-recursive algorithm, parallel computing, graphics processing
unit
https://pp.isofts.kiev.ua
CC BY 4.0
103
Аналітика даних
Вступ
Задача ортогонального (унітарного)
розкладу симетричної матриці до тридіаго-
нальної або стрічкової форми є однією з ба-
зових задач чисельної лінійної алгебри.
Вона виникає у багатьох прикладних зада-
чах, зокрема, при знаходженні власних зна-
чень матриць, розв’язуванні задач спектра-
льного аналізу, моделюванні фізичних про-
цесів та обробці великих масивів даних. У
більшості сучасних алгоритмів обчислення
власних значень попереднє перетворення
матриці до тридіагонального або стрічко-
вого вигляду є необхідним етапом, що зна-
чною мірою визначає загальну обчислюва-
льну складність задачі. Зі зростанням роз-
мірів задач та розвитком високопродуктив-
них обчислювальних систем особливого
значення набуває розробка алгоритмів, які
дозволяють ефективно використовувати
паралельні архітектури. Традиційні методи
розкладу матриць, хоча й забезпечують ви-
соку чисельну стабільність, часто мають
обмежені можливості для паралелізації. У
зв’язку з цим актуальною є задача розробки
нових алгоритмічних підходів, що поєдну-
ють чисельну стійкість класичних методів
із високою ефективністю паралельного ви-
конання. Одним із перспективних напрямів
розв’язання цієї задачі є використання бло-
ково-рекурсивних алгоритмів, які дозволя-
ють природно організувати паралельні об-
числення та ефективно використовувати
сучасні високопродуктивні обчислювальні
системи.
Задача тридіагоналізації матриці ті-
сно пов’язана з класичними методами QR-
розкладу, які широко застосовуються у чи-
сельній лінійній алгебрі завдяки своїй чисе-
льній стійкості. Основні підходи до побу-
дови QR-розкладу — обертання Гівенса,
ортогоналізація Грама–Шмідта та відбиття
Хаусхолдера — детально описані у класич-
них роботах Голуба і Ван Лоана [1], а також
Трефезена і Бау [2]. Одним із важливих на-
прямів розвитку алгоритмів чисельної лі-
нійної алгебри стало використання методів
швидкого множення матриць. У роботі
Штрассена показано можливість зниження
складності множення матриць до 𝑶𝑶(𝒏𝒏𝐥𝐥𝐥𝐥𝐥𝐥𝟐𝟐 𝟕𝟕)
[3], а подальші дослідження Копперсміта і
Вінограда дозволили зменшити асимптоти-
чну оцінку до 𝑶𝑶(𝒏𝒏𝟐𝟐.𝟑𝟑𝟑𝟑) [4]. Ці результати
створили передумови для розробки ефекти-
вних блоково-рекурсивних алгоритмів у за-
дачах чисельної лінійної алгебри.
Однією з ранніх робіт у цьому на-
прямі є дослідження Шонхаге [5], у якому
розглянуто швидке виконання унітарних
перетворень великих матриць. У роботі по-
казано, що такі перетворення можуть бути
організовані таким чином, що їхня асимп-
тотична складність має той самий порядок,
що і складність множення матриць. Пода-
льший розвиток цієї ідеї привів до ство-
рення рекурсивних алгоритмів факториза-
ції матриць. Зокрема, у роботах [6, 7, 8] ро-
зглядалось застосування блоково-рекурси-
вних алгоритмів до задач унітарного роз-
кладу матриць, пов’язаних із задачею QR-
розкладу.
Важливим результатом цього напря-
мку стала робота Деммела та співавторів
[9], у якій доведено, що швидкі алгоритми
множення матриць можуть бути викорис-
тані для розв’язання стандартних задач чи-
сельної лінійної алгебри (SVD, LU, QR) із
тією самою асимптотичною складністю,
зберігаючи водночас чисельну стійкість.
Автори показали, що навіть у разі викорис-
тання так званих слабо стійких алгоритмів
множення матриць забезпечується зворо-
тна стійкість відповідних матричних роз-
кладів.
Подальші дослідження показали, що
рекурсивні матричні алгоритми можуть
ефективно масштабуватися у паралельних
обчислювальних системах із динамічним
розподілом обчислювальних ресурсів [6].
Для реалізації таких алгоритмів розроблено
середовище виконання DAP, яке забезпечує
ефективне розпаралелювання рекурсивних
алгоритмів у системах із розподіленою
пам’яттю [10, 11]. Це створює передумови
для розробки нових блоково-рекурсивних
алгоритмів.
Метою даної роботи є розробка та
дослідження ефективного алгоритму орто-
гонального (унітарного) розкладу симетри-
чних матриць до стрічкової форми на ос-
нові блоково-рекурсивних QR та QP пере-
104
Аналітика даних
творень. Для досягнення цієї мети у роботі
запропоновано алгоритм Q3RP, що базу-
ється на послідовному застосуванні QR та
QP перетворень. Також проведено аналіз
обчислювальної складності алгоритму та
експериментальне дослідження його про-
дуктивності у паралельному середовищі
виконання DAP із використанням приско-
рення матричних операцій на графічному
процесорі.
Ортогональний розклад
симетричної матриці
Ортогональний (унітарний) розклад
симетричної матриці до тридіагональної
форми є одним із ключових етапів у бага-
тьох алгоритмах знаходження власних зна-
чень. Класичні методи такого розкладу, зо-
крема, метод відбиттів Хаусхолдера, мають
складність порядку 𝑶𝑶(𝒏𝒏𝟑𝟑) і реалізуються
переважно як послідовні алгоритми, що об-
межує їхню ефективність у паралельних об-
числювальних системах [1, 2]. У зв’язку з
цим значний інтерес становлять блоково-
рекурсивні алгоритми, які дозволяють
краще використовувати паралельну архіте-
ктуру сучасних обчислювальних платформ.
QR-розклад. Одним із таких підхо-
дів є блоково-рекурсивні алгоритми, що ба-
зуються на QR-розкладі матриці. QR-розк-
лад представляє довільну матрицю 𝑨𝑨 у ви-
гляді добутку ортогональної матриці Q та
верхньої трикутної матриці 𝑹𝑹:
𝑨𝑨 = 𝑸𝑸𝑻𝑻𝑹𝑹, 𝑸𝑸𝑻𝑻𝑸𝑸 = 𝑸𝑸𝑸𝑸𝑻𝑻 = 𝑰𝑰,
(1)
де 𝑸𝑸 — ортогональна матриця, а 𝑹𝑹 — вер-
хня трикутна матриця.
Одним із поширених способів побу-
дови такого розкладу є метод обертань Ґі-
венса, який послідовно занулює піддіагона-
льні елементи матриці за допомогою орто-
гональних матриць повороту. Для двох еле-
ментів вектора (𝒂𝒂, 𝒃𝒃)𝑻𝑻 відповідна матриця
повороту має вигляд
𝒈𝒈(𝛂𝛂, 𝛄𝛄) = ( 𝛂𝛂 𝛄𝛄
−𝛄𝛄 𝛂𝛂) , 𝛂𝛂𝟐𝟐 + 𝛄𝛄𝟐𝟐 = 𝟏𝟏,
(2)
що забезпечує ортогональність перетво-
рення. Узагальненням такого перетворення
для матриці є матриця Ґівенса
𝑮𝑮(𝒊𝒊, 𝒋𝒋) = diag(𝑰𝑰, 𝒈𝒈(𝛂𝛂, 𝛄𝛄), 𝑰𝑰),
(3)
яка відрізняється від одиничної лише у
двох рядках. Завдяки цьому відповідні пе-
ретворення є локальними та придатними
для паралельної реалізації. Для щільної ма-
триці розміру 𝒏𝒏 × 𝒏𝒏 обчислювальна склад-
ність такого алгоритму становить прибли-
зно 𝟐𝟐𝒏𝒏𝟑𝟑 арифметичних операцій.
Подальший розвиток цієї ідеї привів
до появи блоково-рекурсивних алгоритмів
QR-розкладу. У роботі [6] запропоновано
алгоритм, у якому матриця
𝑴𝑴 = (𝑨𝑨 𝑩𝑩
𝑪𝑪 𝑫𝑫)
(4)
розбивається на чотири блоки однакового
розміру. Алгоритм складається з трьох ос-
новних етапів.
На першому етапі виконується QR-
розклад блоку 𝑪𝑪
𝑸𝑸𝟏𝟏𝑪𝑪 = 𝑪𝑪𝑼𝑼,
(5)
у результаті чого формується верхньотри-
кутна матриця 𝑪𝑪𝑼𝑼 та ортогональне перетво-
рення 𝑸𝑸𝟏𝟏. Отримане перетворення застосо-
вується до всієї матриці.
Другим етапом є спеціальна проце-
дура, яка називається QP-розкладом – Рис.
1. Вона використовується для занулення
елементів області, утвореної нижньою три-
кутною частиною блоку 𝑨𝑨 та верхньою три-
кутною частиною блоку 𝑪𝑪𝑼𝑼. Ця область має
форму так званого «паралелограму». У про-
цесі QP-розкладу будується ортогональне
перетворення 𝑸𝑸𝟐𝟐, яке перетворює матрицю
𝑷𝑷 = ( 𝑨𝑨
𝑪𝑪𝑼𝑼)
(6)
до верхньотрикутного вигляду.
У результаті застосування перетво-
рення 𝑸𝑸𝟐𝟐 матриця 𝑴𝑴𝟏𝟏 набуває вигляду
𝑴𝑴𝟐𝟐 = 𝑸𝑸𝟐𝟐 ( 𝑨𝑨 𝑩𝑩
𝑪𝑪𝑼𝑼 𝑫𝑫𝟏𝟏
) = (𝑨𝑨𝑼𝑼 𝑩𝑩𝟏𝟏
𝟎𝟎 𝑫𝑫𝟐𝟐
),
(7)
105
Аналітика даних
де 𝑨𝑨𝑼𝑼 — верхньотрикутний блок, 𝑩𝑩𝟏𝟏 —
оновлений правий блок, а 𝑫𝑫𝟐𝟐 — модифіко-
ваний нижній блок матриці.
На третьому етапі виконується QR-
розклад оновленого блоку 𝑫𝑫𝟐𝟐, що приво-
дить до загального розкладу
𝑴𝑴 = 𝑸𝑸 𝑹𝑹,
(8)
де
𝑸𝑸 = diag(𝑰𝑰, 𝑸𝑸𝟑𝟑) 𝑸𝑸𝟐𝟐 diag(𝑰𝑰, 𝑸𝑸𝟏𝟏).
(9)
QP-розклад. QP-розклад будується
за блоково-рекурсивним принципом. Ми
шукаємо розклад матриці 𝑷𝑷:
𝑷𝑷 = ( 𝑨𝑨
𝑪𝑪𝑼𝑼)
𝑄𝑄2𝑃𝑃 = 𝑃𝑃𝑈𝑈, 𝑃𝑃 = 𝑄𝑄2
𝑇𝑇𝑃𝑃𝑈𝑈, 𝑄𝑄2
−1 = 𝑄𝑄2
𝑇𝑇.
(10)
Матриця 𝑃𝑃 є лівою половиною мат-
риці 𝑀𝑀1. Ми розбили кожен із двох блоків,
що утворюють матрицю 𝑃𝑃 (розміру 2𝑛𝑛 ×
𝑛𝑛), на чотири рівні частини. У результаті
отримано 8 блоків, включаючи один нульо-
вий блок і два верхньотрикутні блоки 𝑓𝑓𝑈𝑈 та
ℎ𝑈𝑈. Необхідно занулити всі елементи мат-
риці P, розташовані між верхньою та ниж-
ньою діагоналями, включно із нижньою ді-
агоналлю. Нас цікавить блокова процедура.
Оскільки 𝑛𝑛 є парним, паралелограм, утво-
рений цими діагоналями, можна розділити
на 4 частини за допомогою двох середніх
ліній. У результаті отримуємо 4 однакові
паралелограми. Для розкладу кожного з
них процедура QP викликається рекурси-
вно 4 рази. Обробка виконується у такому
порядку: спочатку нижній лівий паралелог-
рам (𝑃𝑃𝑙𝑙𝑙𝑙), потім одночасно зануляються вер-
хній лівий (𝑃𝑃𝑢𝑢𝑙𝑙) та нижній правий (𝑃𝑃𝑙𝑙𝑙𝑙), і на
завершення виконується розклад паралело-
грама у верхньому правому куті (𝑃𝑃𝑢𝑢𝑙𝑙). Від-
повідні ортогональні матриці розміру 𝑛𝑛 ×
𝑛𝑛 позначаються 𝑄𝑄𝑙𝑙𝑙𝑙, 𝑄𝑄𝑢𝑢𝑙𝑙, 𝑄𝑄𝑙𝑙𝑙𝑙та𝑄𝑄𝑢𝑢𝑙𝑙.
𝑃𝑃 = (𝐴𝐴𝐶𝐶𝑈𝑈) = (
𝑎𝑎 𝑐𝑐
𝑏𝑏 𝑑𝑑
𝑓𝑓𝑈𝑈 𝑔𝑔
0 ℎ𝑈𝑈
),
𝑃𝑃𝑈𝑈 = (𝐴𝐴𝑈𝑈0) = (
𝑎𝑎𝑈𝑈 𝑐𝑐1
0 𝑑𝑑𝑈𝑈
0 0
0 0
).
𝑃𝑃𝑙𝑙𝑙𝑙 = 𝑄𝑄𝑙𝑙𝑙𝑙 ( 𝑏𝑏 𝑑𝑑
𝑓𝑓𝑈𝑈 𝑔𝑔) = (𝑏𝑏𝑈𝑈 𝑑𝑑1
0 𝑔𝑔1
),
𝑃𝑃𝑢𝑢𝑙𝑙 = 𝑄𝑄𝑢𝑢𝑙𝑙 ( 𝑎𝑎 𝑐𝑐
𝑏𝑏𝑈𝑈 𝑑𝑑1
) = (𝑎𝑎𝑈𝑈 𝑐𝑐1
0 𝑑𝑑2
),
𝑃𝑃𝑙𝑙𝑙𝑙 = 𝑄𝑄𝑙𝑙𝑙𝑙 (0 𝑔𝑔1
0 ℎ𝑈𝑈) = (0 𝑔𝑔𝑈𝑈
0 0 ),
𝑃𝑃𝑢𝑢𝑙𝑙 = 𝑄𝑄𝑢𝑢𝑙𝑙 (0 𝑑𝑑2
0 𝑔𝑔𝑈𝑈) = (0 𝑑𝑑𝑈𝑈
0 0 ).
(11)
У результаті отримуємо матриці 𝑄𝑄2
та 𝑃𝑃𝑈𝑈:
𝑄𝑄2 = 𝑄𝑄𝑢𝑢𝑙𝑙̅̅ ̅̅ ̅ 𝑄𝑄2̅̅ ̅ 𝑄𝑄𝑙𝑙𝑙𝑙̅̅ ̅̅ , 𝑄𝑄2𝑃𝑃 = 𝑃𝑃𝑈𝑈,
(12)
де
𝑄𝑄𝑢𝑢𝑙𝑙̅̅ ̅̅ ̅ = diag (𝐼𝐼𝑛𝑛
2
, 𝑄𝑄𝑢𝑢𝑙𝑙, 𝐼𝐼𝑛𝑛
2
),
𝑄𝑄2̅̅ ̅ = diag(𝑄𝑄𝑢𝑢𝑙𝑙, 𝑄𝑄𝑙𝑙𝑙𝑙),
𝑄𝑄𝑙𝑙𝑙𝑙̅̅ ̅̅ = diag (𝐼𝐼𝑛𝑛
2
, 𝑄𝑄𝑙𝑙𝑙𝑙, 𝐼𝐼𝑛𝑛
2
).
(13)
Схему блоково-рекурсивного QR-
розкладу представлено на Рис. 1.
Рис. 1. Схема блоково-рекурсивного
QR-розкладу
Складність QP-розкладу можна оці-
нити, підрахувавши кількість множень мат-
ричних блоків. У процесі виконання алго-
ритму необхідно виконати 28 множень бло-
ків розміру 𝒏𝒏/𝟐𝟐 × 𝒏𝒏/𝟐𝟐: 8 множень викори-
стовуються у процесі обчислення матриць
𝑷𝑷𝒍𝒍𝒍𝒍 та 𝑷𝑷𝒖𝒖𝒍𝒍, а ще 20 — при побудові матриці
𝑸𝑸𝟐𝟐. Тому загальна кількість операцій задо-
вольняє рекурентному співвідношенню
𝑪𝑪𝒒𝒒𝒒𝒒(𝟐𝟐𝒏𝒏) = 𝟒𝟒𝑪𝑪𝒒𝒒𝒒𝒒(𝒏𝒏) + 𝟐𝟐𝟐𝟐𝑴𝑴(𝒏𝒏/𝟐𝟐),
(14)
де 𝑴𝑴(𝒏𝒏) — складність множення матриць
розміру 𝒏𝒏 × 𝒏𝒏.
106
Аналітика даних
Нехай складність множення двох
матриць розміру 𝒏𝒏 × 𝒏𝒏 дорівнює
𝑴𝑴(𝒏𝒏) = 𝛄𝛄𝒏𝒏𝛃𝛃,
(15)
а 𝒏𝒏 = 𝟐𝟐𝒌𝒌. Тоді рекурентне співвідношення
набуває вигляду
𝑪𝑪𝒒𝒒𝒒𝒒(𝟐𝟐𝒌𝒌+𝟏𝟏) = 𝟒𝟒𝑪𝑪𝒒𝒒𝒒𝒒(𝟐𝟐𝒌𝒌) + 𝟐𝟐𝟐𝟐𝑴𝑴(𝟐𝟐𝒌𝒌−𝟏𝟏).
(16)
Розгортаючи рекурсію, отримуємо
𝑪𝑪𝒒𝒒𝒒𝒒(𝟐𝟐𝒌𝒌+𝟏𝟏) = 𝟒𝟒𝒌𝒌𝑪𝑪𝒒𝒒𝒒𝒒(𝟐𝟐)
+ 𝟐𝟐𝟐𝟐𝛄𝛄 ∑ 𝟒𝟒𝒌𝒌−𝒊𝒊−𝟏𝟏𝟐𝟐𝒊𝒊𝛃𝛃
𝒌𝒌−𝟏𝟏
𝒊𝒊=𝟎𝟎
.
(17)
Для розв'язання рекурентного
співвідношення врахуємо базовий випадок
𝑪𝑪𝒒𝒒𝒒𝒒(𝟐𝟐) = 𝟔𝟔.
(18)
Тоді отримуємо явний вираз для
складності QP-розкладу
𝑪𝑪𝒒𝒒𝒒𝒒(𝒏𝒏) = 𝟐𝟐𝟐𝟐𝛄𝛄(𝒏𝒏𝛃𝛃 − 𝒏𝒏𝟐𝟐)
𝟐𝟐𝛃𝛃(𝟐𝟐𝛃𝛃 − 𝟒𝟒) + 𝟔𝟔𝒏𝒏𝟐𝟐.
(19)
Обмежуючись старшим ступенем за
n, отримуємо асимптотичну оцінку
𝑪𝑪𝒒𝒒𝒒𝒒(𝒏𝒏) = 𝛗𝛗(𝛃𝛃)𝛄𝛄𝒏𝒏𝛃𝛃, 𝛗𝛗(𝛃𝛃) = 𝟐𝟐𝟐𝟐
𝟐𝟐𝛃𝛃(𝟐𝟐𝛃𝛃 − 𝟒𝟒).
𝛗𝛗(𝟑𝟑) = 𝟕𝟕
𝟐𝟐 , 𝛗𝛗(𝐥𝐥𝐥𝐥𝐥𝐥𝟐𝟐 𝟕𝟕) = 𝟒𝟒
𝟑𝟑.
(20)
Складність блоково-рекурсивного
QR-розкладу визначається з урахуванням
того, що алгоритм містить два рекурсивні
виклики QR-процедури, одну процедуру
QP-розкладу та кілька операцій множення
матриць. Якщо позначити складність QP-
розкладу як
𝑪𝑪𝒒𝒒𝒒𝒒(𝒏𝒏) = 𝛂𝛂𝒏𝒏𝛃𝛃,
(21)
а складність множення матриць як (15), то
кількість операцій QR-розкладу задоволь-
няє рекурентному співвідношенню
𝑪𝑪𝒒𝒒𝒒𝒒(𝒏𝒏) = 𝟐𝟐𝑪𝑪𝒒𝒒𝒒𝒒(𝒏𝒏/𝟐𝟐) + 𝑪𝑪𝒒𝒒𝒒𝒒(𝒏𝒏/𝟐𝟐)
+ 𝟐𝟐𝑴𝑴(𝒏𝒏/𝟐𝟐)
(22)
Нехай 𝒏𝒏 = 𝟐𝟐𝒌𝒌. Тоді
𝑪𝑪𝒒𝒒𝒒𝒒(𝟐𝟐𝒌𝒌) = 𝟐𝟐𝑪𝑪𝒒𝒒𝒒𝒒(𝟐𝟐𝒌𝒌−𝟏𝟏) + 𝑪𝑪𝒒𝒒𝒒𝒒(𝟐𝟐𝒌𝒌−𝟏𝟏)
+ 𝟐𝟐𝑴𝑴(𝟐𝟐𝒌𝒌−𝟏𝟏)
= 𝟐𝟐𝒌𝒌𝑪𝑪𝒒𝒒𝒒𝒒(𝟐𝟐𝟎𝟎) + ∑ 𝟐𝟐𝒌𝒌−𝒊𝒊𝑪𝑪𝒒𝒒𝒒𝒒(𝟐𝟐𝒊𝒊−𝟏𝟏)
𝒌𝒌
𝒊𝒊=𝟎𝟎
+ 𝟐𝟐 ∑ 𝟐𝟐𝒌𝒌−𝒊𝒊𝑴𝑴(𝟐𝟐𝒊𝒊−𝟏𝟏)
𝒌𝒌
𝒊𝒊=𝟎𝟎
= 𝛄𝛄(𝛟𝛟 + 𝟐𝟐) ∑ 𝟐𝟐𝒌𝒌−𝒊𝒊(𝟐𝟐𝛃𝛃(𝒊𝒊−𝟏𝟏))
𝒌𝒌
𝒊𝒊=𝟎𝟎
= 𝛄𝛄(𝛟𝛟 + 𝟐𝟐) 𝟐𝟐𝛃𝛃𝒏𝒏𝛃𝛃 − 𝟐𝟐𝒏𝒏
𝟐𝟐𝛃𝛃(𝟐𝟐𝛃𝛃 − 𝟐𝟐)
(23)
Обмежуючись старшим ступенем за
𝒏𝒏, отримуємо:
𝑪𝑪𝒒𝒒𝒒𝒒(𝒏𝒏) = 𝛒𝛒(𝛃𝛃)𝛄𝛄𝒏𝒏𝛃𝛃, 𝛒𝛒(𝛃𝛃)
= (𝛟𝛟 + 𝟐𝟐)/(𝟐𝟐𝛃𝛃 − 𝟐𝟐)
𝛒𝛒(𝟑𝟑) = 𝟕𝟕𝟏𝟏/𝟒𝟒𝟐𝟐, 𝛒𝛒(𝐥𝐥𝐥𝐥𝐥𝐥𝟐𝟐 𝟕𝟕) = 𝟏𝟏𝟒𝟒/𝟗𝟗
(24)
Алгоритм Q3RP. Описані вище вла-
стивості блоково-рекурсивних QR та QP
розкладів дозволяють використовувати їх
як базові операції для побудови більш скла-
дних алгоритмів матричних перетворень.
Одним із них є алгоритм, який ми будемо
називати Q3RP, призначений для ортогона-
льного розкладу симетричної матриці до
стрічкової форми.
Основна ідея алгоритму полягає в
поетапному зануленні цілих блокових ша-
рів симетричної матриці за допомогою пос-
лідовності QR та QP розкладів. На відміну
від класичних алгоритмів, де елементи за-
нулюються послідовними перетвореннями
окремих рядків або стовпців, у цьому алго-
ритмі занулення виконується для цілих бло-
ків матриці. Це дозволяє відокремлювати та
незалежно обчислювати відповідні підза-
дачі, що забезпечує високий рівень парале-
лізму та робить алгоритм ефективним для
обробки великих матриць на кластерах з ро-
зподіленою пам’яттю.
Нехай задано симетричну матрицю
𝑴𝑴 розміру 𝒏𝒏 × 𝒏𝒏. На першому етапі мат-
риця розбивається на чотири великі блоки
107
Аналітика даних
𝑴𝑴 = (𝑨𝑨 𝑩𝑩
𝑪𝑪 𝑫𝑫),
(25)
де кожен блок має розмір 𝒏𝒏/𝟐𝟐 × 𝒏𝒏/𝟐𝟐. Далі
кожен із цих блоків додатково ділиться на
чотири підблоки розміру 𝒏𝒏/𝟒𝟒 × 𝒏𝒏/𝟒𝟒. Та-
ким чином, матриця розбивається на 16 рі-
вних блоків.
𝑴𝑴 = (𝑨𝑨 𝑩𝑩
𝑪𝑪 𝑫𝑫) =
(
𝒂𝒂𝟏𝟏,𝟏𝟏 𝒂𝒂𝟏𝟏,𝟐𝟐 𝒃𝒃𝟏𝟏,𝟏𝟏 𝒃𝒃𝟏𝟏,𝟐𝟐
𝒂𝒂𝟐𝟐,𝟏𝟏 𝒂𝒂𝟐𝟐,𝟐𝟐 𝒃𝒃𝟐𝟐,𝟏𝟏 𝒃𝒃𝟐𝟐,𝟐𝟐
𝒄𝒄𝟏𝟏,𝟏𝟏 𝒄𝒄𝟏𝟏,𝟐𝟐 𝒅𝒅𝟏𝟏,𝟏𝟏 𝒅𝒅𝟏𝟏,𝟐𝟐
𝒄𝒄𝟐𝟐,𝟏𝟏 𝒄𝒄𝟐𝟐,𝟐𝟐 𝒅𝒅𝟐𝟐,𝟏𝟏 𝒅𝒅𝟐𝟐,𝟐𝟐)
(26)
Схематично це представлено на Рис. 2.
Рис. 2. Схема розбиття симетричної мат-
риці на блоки для алгоритму Q3RP
Першим кроком алгоритму є вико-
нання QR-розкладу блоку 𝒄𝒄𝟐𝟐,𝟏𝟏. У результаті
цей блок набуває верхньотрикутної форми,
а отримана ортогональна матриця 𝑸𝑸𝟏𝟏 задає
відповідне ортогональне перетворення, яке
застосовується до інших блоків матриці. За-
вдяки симетричності матриці таке перетво-
рення одночасно впливає і на симетричний
блок 𝒃𝒃𝟏𝟏,𝟐𝟐. Після цього виконується QP-роз-
клад, який дозволяє усунути решту елемен-
тів у відповідному верхньотрикутному
блоці. Добуток отриманих ортогональних
матриць утворює узагальнену матрицю по-
вороту 𝑸𝑸𝟑𝟑. Поточна матриця домножується
на 𝑸𝑸𝟑𝟑𝑻𝑻 зліва та на 𝑸𝑸𝟑𝟑 справа, що забезпечує
збереження її симетричності. На наступ-
ному етапі аналогічно обробляються блоки
𝒄𝒄𝟏𝟏,𝟏𝟏𝑼𝑼 та симетричний до нього блок (𝒃𝒃𝟏𝟏,𝟏𝟏)
𝑼𝑼
за допомогою QP-розкладу. Далі послідо-
вно застосовуються QR- та QP-розклади
для блоків 𝒄𝒄𝟐𝟐,𝟐𝟐′′ та 𝒄𝒄𝟏𝟏,𝟐𝟐′′′ . Кожен із цих кроків
породжує нові ортогональні перетворення,
які комбінуються у відповідну матрицю по-
вороту для поточного шару.
У результаті послідовного застосу-
вання цих перетворень усі блоки матриці 𝑪𝑪
зануляються, що приводить до повного за-
нулення нижнього лівого блоку матриці ра-
зом із відповідними симетричними елемен-
тами верхнього правого блоку.
Завершальним етапом процедури є
QR-розклад блоку 𝒅𝒅𝟐𝟐,𝟏𝟏
(𝟓𝟓). Отримана матриця
повороту 𝑸𝑸𝟖𝟖 знову застосовується до мат-
риці зліва та справа, що завершує форму-
вання стрічкової структури. У результаті
отримуємо представлення початкової мат-
риці у вигляді
𝑺𝑺 = 𝑸𝑸𝑻𝑻𝑴𝑴𝑸𝑸,
(27)
де 𝑸𝑸 = 𝑸𝑸𝟏𝟏𝑸𝑸𝟐𝟐 ⋯𝑸𝑸𝟖𝟖 є добутком ортогональ-
них матриць, отриманих на кожному кроці
алгоритму, а 𝑺𝑺 — результуюча стрічкова
матриця. Матриця 𝑺𝑺 має спеціальну струк-
туру: усі елементи, що знаходяться поза ме-
жами діагональної смуги шириною 𝒑𝒑 = 𝒏𝒏
𝟒𝟒,
дорівнюють нулю. Ширина цієї смуги ви-
значається розміром блоків, на які розбива-
ється початкова матриця. Отримана стріч-
кова форма є проміжним результатом, який
надалі може бути перетворений до тридіа-
гонального вигляду за допомогою спеціалі-
зованих алгоритмів звуження смуги, зок-
рема, алгоритму Ribbon, який буде предста-
влено в наступній статті.
Складність Q3RP. Нехай задано си-
метричну матрицю 𝑴𝑴 розміру 𝟐𝟐𝒌𝒌 × 𝟐𝟐𝒌𝒌. Для
оцінки складності алгоритму використаємо
складності QR та QP розкладів, описані у
попередньому розділі, та припустимо, що
складність множення матриць задається
функцією
𝑴𝑴(𝒏𝒏) = 𝛄𝛄𝒏𝒏𝛃𝛃.
(28)
Аналіз послідовності операцій алго-
ритму Q3RP показує, що для матриці роз-
міру 𝒏𝒏 необхідно виконати:
три QR-розклади для матриць роз-
міру 𝒏𝒏/𝟒𝟒;
108
Аналітика даних
три QP-розклади для матриць роз-
міру 𝒏𝒏/𝟒𝟒;
шість множень матриць розміру 𝒏𝒏/
𝟐𝟐;
двадцять вісім множень матриць ро-
зміру 𝒏𝒏/𝟒𝟒.
Оскільки множення матриць роз-
міру 𝒏𝒏/𝟐𝟐 можна звести до множень блоків
розміру 𝒏𝒏/𝟒𝟒, усі матричні множення можна
представити як множення блоків одного ро-
зміру 𝒏𝒏/𝟒𝟒. Загальна кількість таких мно-
жень становить 50. Тож загальну склад-
ність алгоритму можна записати у вигляді
𝑪𝑪𝒏𝒏 = 𝟑𝟑𝑪𝑪𝒒𝒒𝒒𝒒(𝒏𝒏/𝟒𝟒) + 𝟑𝟑𝑪𝑪𝒒𝒒𝒒𝒒(𝒏𝒏/𝟒𝟒)
+ 𝟓𝟓𝟓𝟓𝟓𝟓(𝒏𝒏/𝟒𝟒)
= (𝟑𝟑𝛟𝛟 + 𝟑𝟑𝛒𝛒 + 𝟓𝟓𝟓𝟓)𝛄𝛄(𝒏𝒏𝛃𝛃)/𝟒𝟒𝜷𝜷
(29)
тобто:
𝑪𝑪𝒏𝒏 = 𝛕𝛕𝛄𝛄𝒏𝒏𝜷𝜷
𝛕𝛕(𝛃𝛃) = (𝟑𝟑𝛟𝛟 + 𝟑𝟑𝛒𝛒 + 𝟓𝟓𝟓𝟓)/𝟒𝟒𝜷𝜷
𝛟𝛟 + 𝛒𝛒 = (𝟐𝟐𝟐𝟐(𝟐𝟐𝛃𝛃 − 𝟏𝟏) + 𝟐𝟐 (𝟐𝟐𝛃𝛃(𝟐𝟐𝛃𝛃 − 𝟒𝟒)))
/ ((𝟐𝟐𝛃𝛃 − 𝟐𝟐) (𝟐𝟐𝛃𝛃(𝟐𝟐𝛃𝛃 − 𝟒𝟒)))
𝛟𝛟 + 𝛒𝛒 = (𝟐𝟐 · 𝟐𝟐𝟐𝟐𝛃𝛃 + 𝟒𝟒 · 𝟐𝟐𝛃𝛃 − 𝟐𝟐𝟐𝟐)
/ ((𝟐𝟐𝛃𝛃 − 𝟐𝟐) (𝟐𝟐𝛃𝛃(𝟐𝟐𝛃𝛃 − 𝟒𝟒)))
𝛕𝛕(𝛃𝛃) = (𝟓𝟓𝟓𝟓 · 𝟐𝟐𝟑𝟑𝛃𝛃 + 𝟒𝟒𝟏𝟏𝟐𝟐 · 𝟐𝟐𝛃𝛃 − 𝟐𝟐𝟒𝟒 − 𝟐𝟐𝟐𝟐𝟐𝟐
· 𝟐𝟐𝟐𝟐𝛃𝛃)
/ ((𝟐𝟐𝛃𝛃 − 𝟐𝟐)𝟐𝟐𝟑𝟑𝛃𝛃(𝟐𝟐𝛃𝛃 − 𝟒𝟒))
𝛕𝛕(𝟑𝟑) = 𝟓𝟓. 𝟗𝟗𝟏𝟏, 𝛕𝛕(𝐥𝐥𝐥𝐥𝐥𝐥𝟐𝟐 𝟐𝟐) = 𝟏𝟏. 𝟐𝟐𝟓𝟓
(30)
Отже, асимптотична складність ал-
горитму Q3RP має той самий порядок, що і
складність множення матриць. Викорис-
тання швидких алгоритмів множення мат-
риць, зокрема, алгоритму Штрассена або
його узагальнень, дозволяє зменшити асим-
птотичний показник степеня 𝛃𝛃 і, відпо-
відно, покращити теоретичну оцінку склад-
ності алгоритму Q3RP. Для порівняння кла-
сичний алгоритм розкладу симетричної ма-
триці до тридіагональної форми на основі
відбиттів Хаусхолдера має складність при-
близно 𝟒𝟒
𝟑𝟑 𝒏𝒏
𝟑𝟑 арифметичних операцій. Однак
цей метод реалізується переважно як послі-
довний алгоритм і не має природної бло-
ково-рекурсивної структури, що обмежує
можливості його ефективної паралельної
реалізації. На відміну від нього, алгоритм
Q3RP використовує блоково-рекурсивні
QR та QP перетворення, що дозволяє вико-
нувати більшість обчислень у вигляді опе-
рацій над матричними блоками. Така стру-
ктура алгоритму добре узгоджується з па-
ралельними обчислювальними середови-
щами та дозволяє ефективно використову-
вати високопродуктивні обчислювальні си-
стеми.
Експерименти
Програмна реалізація алгоритму ви-
конана у середовищі блоково-рекурсивних
паралельних обчислень DAP, яке є децент-
ралізованою системою керування парале-
льними задачами [10, 11]. Система реалізо-
вана мовою Java з використанням MPI для
організації взаємодії між обчислювальними
вузлами. Для запуску алгоритму форму-
ється обчислювальний граф, вершинами
якого є обчислювальні блоки (так звані
дропи), а ребрами — зв’язки передачі даних
між ними. Така структура дозволяє приро-
дно реалізувати рекурсивні матричні алго-
ритми у вигляді набору незалежних підза-
дач, які можуть виконуватися паралельно.
У реалізації алгоритму Q3RP основні обчи-
слювально затратні операції — QR-розк-
лади, QP-розклади та множення матриць —
додатково прискорюються за рахунок вико-
ристання графічного процесора (GPU). Ви-
користання GPU дозволяє ефективно вико-
нувати великі блокові матричні операції,
що є ключовими для блоково-рекурсивної
структури алгоритму. Експерименти про-
водилися на одній обчислювальній ноді з
такими характеристиками: CPU: Intel Core
i5-12600K; RAM: 32 GB DDR4; GPU:
NVIDIA RTX 3070 (8 GB VRAM). Усі обчи-
слення виконувалися у подвійній точності
(double precision). Розмір листової вершини
у середовищі DAP становив 128, що визна-
чає глибину рекурсії алгоритму. Вхідні ма-
триці генерувалися як щільні матриці зі
щільністю 50 % ненульових елементів. Для
кожного розміру матриці час виконання ви-
значався як середнє значення за серією
запусків.
109
Аналітика даних
Таблиця 1.
Час виконання та точність алгоритму
Q3RP
N Час вико-
нання, мс
Абсолютна
похибка
256 98 1.60e−15
512 121 2.34e−15
1024 354 2.71e−15
2048 1863 3.38e−15
4096 9617 4.36e−15
Отримані результати наведені в Таб-
лиці 1 показують, що алгоритм Q3RP де-
монструє стабільну продуктивність при збі-
льшенні розміру матриці.
Рис. 3. Час виконання алгоритму Q3RP
З графіка на Рис. 3 видно, що зі збі-
льшенням розміру задачі час виконання ал-
горитму зростає приблизно кубічно, а це ві-
дповідає теоретичній оцінці складності ал-
горитму.
Крім продуктивності важливою ха-
рактеристикою алгоритму є його чисельна
стабільність. Для оцінки точності було ви-
міряно абсолютну похибку реконструкції
матриці після виконання розкладу. Похи-
бка визначалася як норма різниці між поча-
тковою матрицею та матрицею, відновле-
ною після застосування відповідних орто-
гональних перетворень. Результати показу-
ють, що абсолютна похибка знаходиться на
рівні 10−15, що відповідає машинній точно-
сті для чисел типу double. Це свідчить про
високу чисельну стабільність алгоритму
Q3RP.
Отримані експериментальні резуль-
тати підтверджують ефективність запропо-
нованого алгоритму Q3RP для ортогональ-
ного розкладу симетричних матриць до
стрічкової форми. Блоково-рекурсивна
структура алгоритму дозволяє виконувати
більшість обчислень у вигляді великих ма-
тричних операцій, що добре узгоджується з
архітектурою сучасних обчислювальних
систем.
Використання прискорення на GPU
для виконання QR- та QP-розкладів і опера-
цій множення матриць дозволяє суттєво
зменшити час виконання алгоритму при
зростанні розміру задачі. Водночас алго-
ритм зберігає високу точність обчислень,
що підтверджується малими значеннями
похибки реконструкції.
Висновки
У роботі описано блоково-рекурси-
вні алгоритми ортогонального (унітарного)
розкладу матриць та отримано оцінки їх-
ньої обчислювальної складності. Показано,
що кількість арифметичних операцій визна-
чається алгоритмом множення матриць,
який використовується у відповідних бло-
кових обчисленнях. Водночас кількість
операцій у запропонованих алгоритмах у
будь-якому випадку менша, ніж у двох опе-
раціях множення матриць того самого роз-
міру. Отже, для матриць великого розміру
ефективними є алгоритми, що використо-
вують швидкі методи множення матриць.
Перевагою запропонованих алгори-
тмів у суперкомп’ютерних обчисленнях є
їхня блокова структура. Класичні алгори-
тми, зокрема, алгоритм Хаусхолдера, не до-
зволяють природно відокремлювати неза-
лежні блокові підзадачі, тому розпаралелю-
вання задач великого розміру пов’язане зі
значними труднощами.
Запропоновані блокові алгоритми
усувають ці обмеження. Вони дають змогу
виокремлювати незалежні блокові підза-
дачі та забезпечують паралельний і конвеє-
рний характер обчислювального процесу,
що є особливо важливим для виконання об-
числень на кластерах із розподіленою
пам’яттю.
Література
1. Golub G. H., Van Loan C. F. Matrix Computa-
tions. Baltimore: Johns Hopkins University
Press, 2013. 756 p.
110
Аналітика даних
2. Trefethen L. N., Bau D. Numerical Linear Al-
gebra. Philadelphia: SIAM, 1997.
3. Strassen V. Gaussian elimination is not optimal.
Numerische Mathematik. 1969. Vol. 13. No. 4.
P. 354–356.
4. Coppersmith D., Winograd S. Matrix multipli-
cation via arithmetic progressions. Journal of
Symbolic Computation. 1990. Vol. 9. No. 3. P.
251–280.
5. Schönhage A. Unitäre Transformationen großer
Matrizen. Numerische Mathematik. 1973. Vol.
20. P. 409–417.
6. Malaschonok G. Recursive matrix algorithms,
distributed dynamic control, scaling, stability.
Computer Science and Information Technolo-
gies (CSIT). Yerevan, 2019. P. 112–115.
doi:10.1109/CSITechnol.2019.8895255.
7. Malaschonok G., Ivaskevich A. Quick recursive
QR decomposition. Proceedings of the Confer-
ence on Mathematical Foundations of Informat-
ics (MFOI-2020). Kyiv, 2020. P. 1–8.
8. Першута П. В. Алгоритм зведення
симетричної матриці до тридіагональної
форми із застосуванням QR та QP-розкладів:
магістерська робота. Київ, 2024. 55 с.
9. Demmel J., Dumitriu I., Holtz O., Kleinberg R.
Fast Linear Algebra is Stable. Numerische
Mathematik. 2007. Vol. 108. P. 59–91.
10.Сідько А. А. Середовище виконання для
блоково-рекурсивних матричних алгоритмів
на суперкомп’ютері з розподіленою
пам’яттю: дис. … д-ра філософії. Київ, 2024.
11.Малашонок Г. І., Сідько А. А. Розподілені
обчислення: ДАП-технологія розпарале-
лювання рекурсивних алгоритмів. Наукові
записки НаУКМА. 2018. № 1. С. 25–32.
References
1. Golub, G.H. and Van Loan, C.F., 2013. Matrix
Computations. Baltimore: Johns Hopkins Uni-
versity Press.
2. Trefethen, L.N. and Bau, D., 1997. Numerical
Linear Algebra. Philadelphia: SIAM.
3. Strassen, V., 1969. Gaussian elimination is not
optimal. Numerische Mathematik, 13(4),
pp.354–356.
4. Coppersmith, D. and Winograd, S., 1990. Ma-
trix multiplication via arithmetic progressions.
Journal of Symbolic Computation, 9(3),
pp.251–280.
5. Schönhage, A., 1973. Unitäre Transformatio-
nen großer Matrizen. Numerische Mathematik,
20, pp.409–417.
6. Malaschonok, G., 2019. Recursive matrix algo-
rithms, distributed dynamic control, scaling,
stability. In: Computer Science and Information
Technologies (CSIT). Yerevan, pp.112–115.
https://doi.org/10.1109/CSITechnol.2019.8895
255
7. Malaschonok, G. and Ivaskevich, A., 2020.
Quick recursive QR decomposition. In: Pro-
ceedings of the Conference on Mathematical
Foundations of Informatics (MFOI-2020).
Kyiv, pp.1–8.
8. Pershuta, P.V., 2024. Algorithm for reducing a
symmetric matrix to tridiagonal form using QR
and QP decompositions (Master’s thesis). Kyiv
(in Ukrainian).
9. Demmel, J., Dumitriu, I., Holtz, O. and Klein-
berg, R., 2007. Fast linear algebra is stable.
Numerische Mathematik, 108, pp.59–91.
10.Sidko, A.A., 2024. Execution environment for
block-recursive matrix algorithms on a distrib-
uted-memory supercomputer (PhD disserta-
tion). Kyiv (in Ukrainian).
11.Malashonok, H.I. and Sidko, A.A., 2018. Dis-
tributed computing: DAP-technology for paral-
lelization of recursive algorithms. Naukovi
zapysky NaUKMA, 1, pp.25–32 (in Ukrainian).
Дата першого надходження до видання:
13.03.2026
Внутрішня рецензія отримана:
16.03.2026
Зовнішня рецензія отримана: 17.03.2026
Дата прийняття статті до друку:
19.03.2026
Дата публікації: 16.04.2026
Про авторів:
1Малашонок Геннадій Іванович,
Доктор фізико-математичних наук,
професор
1 Malaschonok Gennadiy,
Ph.D (doctor, physical and mathematical
sciences), professor
https://orcid.org/0000-0002-9698-6374
2Сухарський Сергій Сергійович,
PhD студент
2 Sukharskyi Sergiy,
Post-graduate student
https://orcid.org/0000-0002-5873-984X
111
Аналітика даних
Місце роботи авторів:
1Національний університет
«Києво-Могилянська академія»
1 National University
“Kyiv–Mohyla Academy”
malaschonok@gmail.com
2Інститут Програмних Систем
НАН України
2 Institute of Software Systems.
National Academy of Sciences of Ukraine
serhii.sukharskyi@gmail.com
|
| id | pp_isofts_kiev_ua-article-895 |
| institution | Problems in programming |
| keywords_txt_mv | keywords |
| language | Ukrainian |
| last_indexed | 2026-04-24T01:00:16Z |
| publishDate | 2026 |
| publisher | PROBLEMS IN PROGRAMMING |
| record_format | ojs |
| resource_txt_mv | ppisoftskievua/ba/eeed49dea8f6f011905f73a42e6cdbba.pdf |
| spelling | pp_isofts_kiev_ua-article-8952026-04-23T22:26:13Z A block-recursive approach to unitary matrix decomposition Блоково-рекурсивний підхід до унітарного розкладу матриць Malaschonok, G.I. Sukharskyi, S.S. unitary matrix decomposition; block-recursive algorithm; parallel computing; graphics processing unit UDC 519.61:004.272 унітарний розклад матриці; блоково-рекурсивний алгоритм; паралельні обчислення; графічний процесор УДК 519.61:004.272 Singular value decomposition is one of the fundamental tools of numerical linear algebra and is widely used in data processing, optimization, computer vision, and machine learning. The classical approach to its computation is based on the Householder algorithm, which reduces an arbitrary matrix to a tridiagonal or bidiagonal form followed by the computation of singular values. However, as the matrix size increases, significant difficulties arise in efficiently parallelizing such algorithms, which limits their performance on modern high-performance computing systems. One promising approach to overcoming these limitations is the use of block-recursive algorithms. Such algo rithms allow the original problem to be decomposed into independent block subproblems that can be processed in parallel, providing a high degree of scalability on multiprocessor systems and distributed-memory clusters. This paper proposes a new block-recursive algorithm for the orthogonal (or unitary) decomposition of symmetric matrices based on QR and QP decompositions. The algorithm constructs orthogonal (unitary) factors and pro duces a band matrix, reducing the computational complexity of subsequent stages of spectral decomposition. A complexity analysis is presented showing that the asymptotic order of the algorithm is determined by the complexity of matrix multiplication.Experimental performance studies were carried out in a parallel computing environment using GPU acceleration for matrix operations. The obtained results demonstrate the efficiency of the proposed algorithm as the problem size increases and confirm the agreement between the experimental com plexity and the theoretical estimates. The reconstruction error of the matrix remains at the level of machine precision for double-precision floating-point numbers, which indicates the numerical stability of the proposed approach.Problems in programming 2026; 1: 102-111 Сингулярне розкладання матриці є одним з основних інструментів чисельної лінійної алгебри та широко застосовується у задачах обробки даних, оптимізації, комп’ютерного зору та машинного навчання. Кла сичний підхід до його обчислення базується на алгоритмі Хаусхолдера, який дозволяє звести довільну матрицю до тридіагональної або бідіагональної форми з подальшим обчисленням сингулярних значень. Однак зі зростанням розмірності матриць виникають суттєві труднощі з ефективним розпаралелюванням таких алгоритмів, що обмежує їхню продуктивність на сучасних високопродуктивних обчислювальних системах. Одним із перспективних підходів до подолання цих обмежень є використання блоково-рекурсивних ал горитмів. Такі алгоритми дозволяють розділяти початкову задачу на незалежні блокові підзадачі, що мо жуть виконуватися паралельно, забезпечуючи високий ступінь масштабованості на багатопроцесорних системах та кластерах із розподіленою пам’яттю. У даній роботі запропоновано новий алгоритм ортогонального (унітарного) блоково-рекурсивного роз кладу симетричних матриць, який базується на використанні QR- та QP-розкладів. Алгоритм дозволяє побудувати ортогональні (унітарні) множники та отримати стрічкову матрицю, що зменшує обчислюва льну складність подальших етапів спектрального розкладу. Отримано оцінку обчислювальної складності алгоритму та показано, що її асимптотичний порядок ви значається складністю множення матриць. Експериментальні дослідження продуктивності проведено у середовищі паралельних обчислень із використанням прискорення матричних операцій на графічному процесорі. Отримані результати демонструють ефективність запропонованого алгоритму при зростанні розміру задачі та підтверджують узгодженість експериментальної складності з теоретичними оцінками. Значення похибки реконструкції матриці знаходяться на рівні машинної точності для чисел подвійної точності, що свідчить про чисельну стабільність запропонованого підходу.Problems in programming 2026; 1: 102-111 PROBLEMS IN PROGRAMMING ПРОБЛЕМЫ ПРОГРАММИРОВАНИЯ ПРОБЛЕМИ ПРОГРАМУВАННЯ 2026-04-23 Article Article application/pdf https://pp.isofts.kiev.ua/index.php/ojs1/article/view/895 PROBLEMS IN PROGRAMMING; No 1 (2026); 102-111 ПРОБЛЕМЫ ПРОГРАММИРОВАНИЯ; No 1 (2026); 102-111 ПРОБЛЕМИ ПРОГРАМУВАННЯ; No 1 (2026); 102-111 1727-4907 uk https://pp.isofts.kiev.ua/index.php/ojs1/article/view/895/949 Copyright (c) 2026 PROBLEMS IN PROGRAMMING |
| spellingShingle | unitary matrix decomposition block-recursive algorithm parallel computing graphics processing unit UDC 519.61:004.272 Malaschonok, G.I. Sukharskyi, S.S. A block-recursive approach to unitary matrix decomposition |
| title | A block-recursive approach to unitary matrix decomposition |
| title_alt | Блоково-рекурсивний підхід до унітарного розкладу матриць |
| title_full | A block-recursive approach to unitary matrix decomposition |
| title_fullStr | A block-recursive approach to unitary matrix decomposition |
| title_full_unstemmed | A block-recursive approach to unitary matrix decomposition |
| title_short | A block-recursive approach to unitary matrix decomposition |
| title_sort | block-recursive approach to unitary matrix decomposition |
| topic | unitary matrix decomposition block-recursive algorithm parallel computing graphics processing unit UDC 519.61:004.272 |
| topic_facet | unitary matrix decomposition block-recursive algorithm parallel computing graphics processing unit UDC 519.61:004.272 унітарний розклад матриці блоково-рекурсивний алгоритм паралельні обчислення графічний процесор УДК 519.61:004.272 |
| url | https://pp.isofts.kiev.ua/index.php/ojs1/article/view/895 |
| work_keys_str_mv | AT malaschonokgi ablockrecursiveapproachtounitarymatrixdecomposition AT sukharskyiss ablockrecursiveapproachtounitarymatrixdecomposition AT malaschonokgi blokovorekursivnijpídhíddounítarnogorozkladumatricʹ AT sukharskyiss blokovorekursivnijpídhíddounítarnogorozkladumatricʹ AT malaschonokgi blockrecursiveapproachtounitarymatrixdecomposition AT sukharskyiss blockrecursiveapproachtounitarymatrixdecomposition |