Parallel non-negative sparse extra-large matrix factorization
This paper proposes parallel methods of non-negative huge sparse matrix factorization. The implementation of proposed methods was tested and compared on huge matrix.Prombles in programming 2014; 2-3: 66-71
Saved in:
| Date: | 2025 |
|---|---|
| Main Author: | |
| Format: | Article |
| Language: | Ukrainian |
| Published: |
PROBLEMS IN PROGRAMMING
2025
|
| Subjects: | |
| Online Access: | https://pp.isofts.kiev.ua/index.php/ojs1/article/view/695 |
| Tags: |
Add Tag
No Tags, Be the first to tag this record!
|
| Journal Title: | Problems in programming |
| Download file: | |
Institution
Problems in programming| id |
pp_isofts_kiev_ua-article-695 |
|---|---|
| record_format |
ojs |
| resource_txt_mv |
ppisoftskievua/04/43a2770240e37002af11fc32f680ca04.pdf |
| spelling |
pp_isofts_kiev_ua-article-6952025-04-09T22:22:32Z Parallel non-negative sparse extra-large matrix factorization Паралелізація невід’ємної факторизації розріджених матриць надвеликої розмірності Nasirov, E.M. UDC 681.3 УДК 681.3 This paper proposes parallel methods of non-negative huge sparse matrix factorization. The implementation of proposed methods was tested and compared on huge matrix.Prombles in programming 2014; 2-3: 66-71 У роботі описано побудову моделі паралелізації обчислення невід’ємної факторизації розріджених матриць надвеликої розмірності. Реалізації запропонованих моделей були порівняні в обробці надвеликої матриці.Prombles in programming 2014; 2-3: 66-71 PROBLEMS IN PROGRAMMING ПРОБЛЕМЫ ПРОГРАММИРОВАНИЯ ПРОБЛЕМИ ПРОГРАМУВАННЯ 2025-04-09 Article Article application/pdf https://pp.isofts.kiev.ua/index.php/ojs1/article/view/695 PROBLEMS IN PROGRAMMING; No 2-3 (2014); 66-71 ПРОБЛЕМЫ ПРОГРАММИРОВАНИЯ; No 2-3 (2014); 66-71 ПРОБЛЕМИ ПРОГРАМУВАННЯ; No 2-3 (2014); 66-71 1727-4907 uk https://pp.isofts.kiev.ua/index.php/ojs1/article/view/695/747 Copyright (c) 2025 PROBLEMS IN PROGRAMMING |
| institution |
Problems in programming |
| baseUrl_str |
https://pp.isofts.kiev.ua/index.php/ojs1/oai |
| datestamp_date |
2025-04-09T22:22:32Z |
| collection |
OJS |
| language |
Ukrainian |
| topic |
UDC 681.3 |
| spellingShingle |
UDC 681.3 Nasirov, E.M. Parallel non-negative sparse extra-large matrix factorization |
| topic_facet |
UDC 681.3 УДК 681.3 |
| format |
Article |
| author |
Nasirov, E.M. |
| author_facet |
Nasirov, E.M. |
| author_sort |
Nasirov, E.M. |
| title |
Parallel non-negative sparse extra-large matrix factorization |
| title_short |
Parallel non-negative sparse extra-large matrix factorization |
| title_full |
Parallel non-negative sparse extra-large matrix factorization |
| title_fullStr |
Parallel non-negative sparse extra-large matrix factorization |
| title_full_unstemmed |
Parallel non-negative sparse extra-large matrix factorization |
| title_sort |
parallel non-negative sparse extra-large matrix factorization |
| title_alt |
Паралелізація невід’ємної факторизації розріджених матриць надвеликої розмірності |
| description |
This paper proposes parallel methods of non-negative huge sparse matrix factorization. The implementation of proposed methods was tested and compared on huge matrix.Prombles in programming 2014; 2-3: 66-71 |
| publisher |
PROBLEMS IN PROGRAMMING |
| publishDate |
2025 |
| url |
https://pp.isofts.kiev.ua/index.php/ojs1/article/view/695 |
| work_keys_str_mv |
AT nasirovem parallelnonnegativesparseextralargematrixfactorization AT nasirovem paralelízacíânevídêmnoífaktorizacíírozrídženihmatricʹnadvelikoírozmírností |
| first_indexed |
2025-07-17T09:53:15Z |
| last_indexed |
2025-07-17T09:53:15Z |
| _version_ |
1850410116272095232 |
| fulltext |
Паралельне програмування. Розподілені системи і мережі
© Е.М. Насиров, 2014
66 ISSN 1727-4907. Проблеми програмування. 2014. № 2–3. Спеціальний випуск
УДК 681.3.
ПАРАЛЕЛІЗАЦІЯ НЕВІД’ЄМНОЇ ФАКТОРИЗАЦІЇ
РОЗРІДЖЕНИХ МАТРИЦЬ НАДВЕЛИКОЇ РОЗМІРНОСТІ
Е.М. Насиров
Київський національний університет імені Тараса Шевченка,
03680, Київ, проспект Академіка Глушкова 4д.
E-mail: enasirov@gmail.com
У роботі описано побудову моделі паралелізації обчислення невід’ємної факторизації розріджених матриць надвеликої розмірності.
Реалізації запропонованих моделей були порівняні в обробці надвеликої матриці.
This paper proposes parallel methods of non-negative huge sparse matrix factorization. The implementation of proposed methods was tested
and compared on huge matrix.
Сьогодні невід’ємна факторизація матриць та тензорів є дуже популярною технологією в штучному інте-
лекті взагалі, та зокрема в комп’ютерній лінгвістиці. Використовуючи невід’ємну факторизацію в рамках пара-
дигми латентно-семантичного аналізу, комп’ютерні лінгвісти застосовують даний підхід для розв’язання таких
прикладних задач, як класифікація, кластеризація текстів та термінів [1, 2], побудова мір семантичної близькос-
ті [3], автоматичне виділення лінгвістичних структур та відношень (Selectional Preferences [4] and Verb Sub-
Categorization Frames [5]) та багато інших.
Дана робота описує побудову моделі паралелізації обчислення невід’ємної факторизації розріджених ма-
триць надвеликої розмірності, що представляє особливий інтерес та актуальність для її застосування у великих
NLP системах загального тематичного напрямку, не обмежених використанням лише для вузьких предметних
областей.
Задача невід’ємної факторизації розріджених матриць надвеликої розмірності постала в процесі розробки
системи визначення міри семантичної близькості-зв’язності за технологією латентного семантичного аналізу
[6]. Для побудови системи був оброблений великий текстовий корпус статей англомовної Вікіпедії. Обробка
корпусу полягала у лексичному аналізі та лематизації лексем речень статей та в обчисленні частот вживання
множини слів та словосполучень англійської мови в складі різних статей англомовної Вікіпедії. В результаті
була побудована велика матриця Слова×Статті яка містила частотну оцінку вживання слів у текстах статей Ві-
кіпедії. Розмірність матриці дорівнює 2,437,234 слів-словосполучень на 4,475,180 статей англомовної Вікіпедії.
Після цього для усунення з матриці випадкових даних був встановлений пороговий рівень Т=3, частотні оцінки
в матриці, які менші за пороговий рівень Т, були обнулені). В результаті була побудована розріджена матриця
надвеликого розміру – 156,236,043 ненульових елементів при розмірі 2,437,234×4,475,180. Для невідємної фак-
торизації розрідженої матриці такого розміру знадобилася розробка спеціальної моделі паралелізації матричних
обчислень, яка була реалізована із застосуванням паралельних розподілених обчислень та обчислень на GPU.
Низка реалізацій факторизації матриць була вже реалізова і представлена у [7 – 10]. Але жодна з розроблених
моделей не є допустимим для поставленої задачі. Деякі з них [7–9] не задовольняють потреби розмірності мат-
риці. Модель запропонована в [10] проводить факторизацію матриць запропонованого розміру за задовільний
час, але потребує занадто великих розрахункових комп’ютерних потужностей, що є не завжди доступним.
Алгоритм невід’ємної матричної факторизації
Алгоритм невід’ємної матричної факторизації [11] виконує декомпозицію невід’ємної матриці V на не-
від’ємні матриці W та H таким чином, щоб:
.HWV
Як функція оцінки може бути використана функція виміру відстані між двома невід’ємними матриця-
ми. Однією з таких мір є квадрат Евклідової метрики:
ji
jiji BABA .
22
Така цільова функція обмежена знизу. Нижня границя 0 досягається тоді і тільки тоді коли
.BA
Паралельне програмування. Розподілені системи і мережі
67
Отже, при використанні Евклідової метрики, факторизація матриці полягає у мінімізації
2
HWV при
умові невід’ємності W та H .
Така цільова функція не зростаюча при наступних правилах:
.
)(
)(
ji
T
ji
T
jiji
WHW
VW
HH (1)
.
)(
)(
ji
T
ji
T
jiji
HWH
VH
WW (2)
Виконання ітерацій алгоритму продовжується дот тих пір, поки не буде досягнута стаціонарна точка, або
не буде виконана максимальна кількість ітерацій.
Аналіз моделі
В табл. 1 представлено необхідні обсяги пам’яті для зберігання матриць W та H при різних k при вико-
ристанні типу даних float (32bit).
Таблиця 1. Необхідні обсяги пам’яті для зберігання матриць W та H при різних k
K 100 200 300
W 0.98Gb 1.95Gb 2.92Gb
H 1.79Gb 3.58Gb 5.37Gb
Всього 2.76Gb 5.53Gb 8.29Gb
Описаний алгоритм на кожній ітерації потребує в 2 рази більше пам’яті для збережень матриць (не вра-
ховуючи потреби на збереження початкової матриці V). Враховуючи такі потреби в пам'яті, для виконання ал-
горитму на одному комп'ютері необхідно проводити збереження даних на жорсткий диск. Далі буде розглянуто
та порівняно 2 підходи до паралелізації алгоритму – це локальний та розподілений.
Паралелізація алгоритму з використанням GPU
Для спрощення, зробимо заміну в (1) та (2): THH . Отримаємо
jit
T
tt
jit
T
jitjit
WWH
WV
HH
)(
)(
)()(
111
1
1
, (4)
jit
T
tt
jit
jitjit
HHW
HV
WW
)(
)(
)()(
1
1
. (5)
Таким чином, обидві формули (4) та (5) ми можемо представити в одному вигляді, завдяки заміні H , W
та TV , або W, H та V на A, B та S відповідно:
.
)(
)(
ij
T
ij
ijij
BBA
BS
AA (6)
Формула (6) може бути представлена як послідовність чотирьох кроків (7):
,BSC (7’)
,BBK T (7”)
,KAD (7”’)
.
ij
ij
ijij
D
C
AA (7””)
Такий порядок обчислень є оптимальним для виразу (6). Ці кроки маюсь обчислювальну складність
)(),()),)((( 22 nkOmkOnSnnzkO та )( nkO відповідно, де )(Snnz – кількість ненульових елементів матри-
Паралельне програмування. Розподілені системи і мережі
68
ці S. Перші три кроки(7’–7”’) підтримуються бібліотеками CUDA cuSPARSE [12] та cuBLAS [13]. Четвертий
крок (7””) для виконання потребує реалізації власного ядра GPU, але в той же час, це відносно швидка операція
і тому може бути виконана на CPU.
Так як матриці занадто великі для збереження в пам’яті GPU, ці операції (7’–7””) мають виконуватись
над частинами, з врахуванням необхідності зменшення обсягів обміну даними з графічним адаптером.
Для виразу (7’) мають розкласти матриці T
tSSSS )( 21 та )( 21 rBBBB , та підрахувати С:
.
1
111
rtt
r
BSBS
BSBS
C
(8)
Для виразу (7”) варто розглядати )( 21 tBBBB , а отже )( 11 t
T
t
T
BBBBK . Підрахунок цього ви-
разу не потребує додаткового завантаження будь-яких матриць у пам'ять графічного адаптера після виконання
(7’).
Для підрахунку (7’”) ми можемо залишити в пам'яті GPU матрицю К і помножити матрицю А як стовп-
чик блоків.
Також можна зменшити використання пам’яті завдяки комбінуванню (7””) та (7”’), так як (7””) – по еле-
ментне правило і єдине інше правило, в якому використовується матриця А – (7”’). Не буде потреби зберігати в
пам’яті матрицю D, якщо виконувати (7””) на блоці з матриці А, який необхідний для обчислення певного бло-
ку матриці D.
Складність підрахунку Евклідової метрики між початковою та факторизовано моделлю складає
)( mnO . Такий підрахунок просто реалізовується з використанням графічних адаптерів.
Розподілений алгоритм
Наступним кроком для покращення швидкодії є використання мережі ПК для проведення обчислень. Ро-
зглянемо можливі розподіленні моделі. Припустимо, що мережа складається з двох вузлів. Можливі наступні
моделі розподілення обчислень.
1. Підраховувати W та H окремо на різних вузлах. Таким чином обидва вузли будуть перебувати в од-
ному з двох альтернативних режимів. Або будуть вузлом підтримки іншого вузла (передача даних на інший
вузол) або будуть проводити підрахунки (обчислювати матрицю за яку відповідають, за допомогою даних,
отриманих з вузла підтримки). В такій моделі розподілення на кожній ітерації ми повинні будемо передати по
мережі дві матриці такого ж розміру, як і W та H. Також вузол розрахунків буде в основному простоювати, тому
що (7а) є найбільш витратним кроком з усіх чотирьох.
2. Розділити матриці W та H на блоки та розподілити їх між вузлами. )( 21 HHH та )( 21 WWW .
Робимо перший вузол відповідальним за 11,WH , другий – за 22 ,WH . При такому підході кожен вузол виступає
як у ролі вузла підтримки, так і у ролі вузла розрахунків одночасно. В такій моделі вузли мають передавати по
мережі обсяг даних рівний ))()((5,1 HofsizeWofsize , де )(Xofsize – обсяг пам’яті, необхідний для зберіган-
ня матриці Х
3. Розділити матриці W та H на блоки та розподілити їх між вузлами. THHH )( 21 та
TWWW )( 21 . Робимо перший вузол відповідальним за 11, WH , другий – за 22 , WH . В цій моделі, аналогічно
до попередньої, кожен вузол мережі виступає одночасно як у ролі вузла підтримки, так і у ролі вузла розрахун-
ків. Обсяг передачі даних вузлом рівний )()( HofsizeWofsize .
В кожній з представлених моделей також є необхідність передачі одної або декількох матриць ];[ kk , але
їх розміром можна знехтувати в порівнянні з матрицями W та H. Для підрахунку функції оцінки в першій та
третій моделі вузлам необхідно передати дані в обсязі
2
))()(( KHofsizeWofsize
,
де K – кількість вузлів в мережі, у другій моделі необхідно передати в )1( K раз більше даних.
Варто також зазначити, що ріст мережі призведе до експоненціального росту сумарного обсягу переда-
них даних, хоча обсяг передачі даних одного вузла буде не більш ніж ))()((2 HofsizeWofsize .
Для кращого розподілення розрахунків між вузлами мережі, враховуючи розрідженість початкової мат-
риці V, необхідно виконати перестановки рядків та стовпчиків для нормалізації кількості ненульових елементів
у кожному блоці.
Паралельне програмування. Розподілені системи і мережі
69
Очевидно, що третя модель найкраще підходить з точки зору мережевого обміну даними та використан-
ня GPU. Тому саме ця модель була використана в нашому дослідженні, реалізована з використання мережі з
чотирьох вузлів і описана далі.
Згідно з моделлю 3 необхідно розділити матриці HW , і V наступним чином:
TWWWWW )( 4321 ,
THHHHH )( 4321 , .
4441
1411
VV
VV
V
Алгоритм складається з трьох основних фаз: ініціалізація, ітерації, підрахунок метрики. Під час фази іні-
ціалізації матриці W, H і V розподіляються між чотирма вузлами так, щоб i-ий вузол отримав матриці ii HW ,
та ,, ikki VV де 4,...,1k . Дана фаза показана на рис. 1, а.
Фаза ітерації складається з двох аналогічних етапів, один для підрахунку H , а другий для W. Кожен
етап включає 3 кроки.
На першому кроці кожний вузол підраховує матрицю T
ii WW розміру kk та надсилає її до агрегато-
ра. Агрегатор об’єднує всі отримані блоки в одну матрицю wK та надсилає отриманий результат всім вузлам.
Цей крок показано на рис. 1, а.
На другому кроці кожен вузол підраховує власне i
TT
i
T
i
T
i
T
i WVVVV )( 4321 . Отримана матриця збігається за
розмірами з H . Далі кожний вузол ділить цю матрицю згідно з початковим розбиттям матриці H та передає
отримані блоки відповідним вузлам. Цей крок показано на рис. 1, b.
На третьому кроці вузли підраховують матрицю wi KH та виконує оновлення матриці H . Цей крок не
потребує передачі даних по мережі.
Ці три кроки необхідні для підрахування матриці H . Після поновлення H , аналогічні кроки необхідно
виконати для W. А саме необхідно підрахувати наступні матриці: T
ii HH , i
TT
i
T
i
T
i
T
i HVVVV )( 4321 та Hi KW .
У фазі підрахунку метрики кожен вузол передає відповідний блок H всім іншим блокам. Після отри-
мання всіх блоків, кожний вузол підраховує відповідну части метрики. Ця фаза також показана на рис. 1, b.
a b
Рис. 1. Розподілена модель з чотирма вузлами
Результати аналізу
Був реалізований розподілений алгоритм з використанням GPU, який був описаний вище, та проведена
факторзація розрідженої матриці оцінок части вживання слів англійської мови в різних статтях Вікіпедії. Для
порівняння швидкодії була також реалізована локальна модель з використанням GPU та жорсткого диску для
запису даних.
Обидві моделі працювали з однаковими вхідними даними. Локальна версія була запущена на одному з
вузлів з тими ж обмеженнями пам'яті. Для проведення тестів були використані наступні апаратні засоби: Intel
Core i7 CPU, NVIDIA GeForce GTX560 1Gb, 8Gb RAM, 1Gbit LAN, SATA III HD.
В табл. 2 порівняно час та ресурси необхідні для виконання однієї ітерації у локальній та розподіленій
версії. В табл. 3 порівняно час та ресурси, необхідні для підрахунку функції оцінки . Дані таблиць отриманні
при 300k .
Паралельне програмування. Розподілені системи і мережі
70
Таблиця 2. Результати порівняння виконання
ітерації локальної та розподіленої версії
Local Distributed
Прочитано 34.44Gb 6.22Gb
Записано 16.58Gb 6.22Gb
Час ітерації
(розрахунки)
58s 62s
Час ітерації (I/O) 729s 287s
Таблиця 3. Результати порівняння обчислення
функції оцінки
Local Distributed
Прочитано 13.66Gb 6.22Gb
Записано 0 6.22Gb
Час (розрахунки) 45865s 11371s
Час (I/O) 192s 280s
На рис. 2 показано графік отриманої залежність як результати факторизації від k та кількості ітерацій.
Рис. 2. Залежність збіжності алгоритму факторизації від k та кількості ітерацій
Застосування моделі до факторизації тензорів
Для того, щоб обчислити невід'ємні матриці-компоненти {A, B, C} зазвичай застосовується обмежений
оптимізаційний підхід для мінімізації підходящої функції оцінки. Зазвичай мінімізують наступну функцію:
2 2 2 2
( [ , , ]) [ , , ] ,F A B CF F F F
D Y A B C Y A B C A B C
де
A ,
B ,
C – невід’ємні регуляційні параметри.
Існує щонайменше три різні підходи до розв’язання такої оптимізаційної задачі. Перший підхід полягає у
використанні векторної форми функції оцінки: 0],,[ CBAYvecXJ . Для розв’язання використовується
метод найменших квадратів. Цей метод для факторизації тензорів вперше використав Паатеро. Якобіан такої
функції може мати розмір QTIITQJ , а отже такий підхід потребує значних розрахункових витрат.
У другому підході Акар, Колда та Дунлаві запропонували штучно оптимізувати функцію оцінки викори-
стовуючи техніку нелінійної зв’язної градієнтної оптимізації [1].
Третім, і найпоширенішим підходом є використання техніки ALS. В цьому підході підраховується граді-
єнт функції оцінки для кожної матриці окремо.
,
,
.
Прирівнюючи кожну компоненту градієнта до 0 та накладаючи умову невід’ємності можна отримати
ефективні та прості правила оновлення матриць:
,
,
.
Паралельне програмування. Розподілені системи і мережі
71
Основними перевагами такого підходу є висока швидкість збіжності та можливість розподілення для за-
дач великої розмірності.
Запропонована вище розподілена модель може бути застосована для факторизації тензорів. Правила оно-
влення є однаковими. Тому буде показано підрахунок на прикладі однієї матриці:
.
Для мережі з двох вузлів необхідно розділити матриці W та H на блоки таким чином, що
та , та розподілити їх між вузлами так, щоб перший вузол був відповідальним за , другий за
.
При такому розподіленні на першому кроці кожний вузол підраховує матрицю розміру kk та
надсилає її до агрегатора. Агрегатор об’єднує всі отримані блоки в одну матрицю CK та надсилає відповідні
частини до вузлів. Далі аналогічним чином виконуються розрахунки для )( BB . На наступному кроці вузлам
необхідно підрахувати та передати відповідні матриці на агрегаті на якому під-
рахується 1S . Після чого вузлам необхідно підрахувати відповідні частини результуючої матирці у вигляді
.
Таким чином можливе використання запропонованого розподіленої моделі до факторизації тензорів.
Висновки
В роботі описано локальний та розподілений підходи до невід’ємної факторизації надвеликих матиць. В
результаті роботи було отримано високі показники швидкодії роботи розподіленого алгоритму в порівнянні з
локальним. Експерименти показали, що процес підрахунку матриць збігається за 100 ітерацій, а час необхідний
для фактризації опрацьованої матриці для розподіленої моделі – 9,6 годин та 21 година для локальної моделі з
використанням GPU. Подальша робота полягає в аналізі та зменшенні обсягів передачі даних між вузлами, а
отже, і часу що для цього необхідний.
1. Xu, W., Liu, X., & Gong, Y. (). Document-clustering based on 71n-negative matrix factorization // In Proceedings of SIGIR’03, July 28–August
1, Toronto – 2003. – P. 267–273.
2. Farial Shahnaz, Michael W. Berry, V. Paul Pauca, Robert J. Plemmons Document clustering using nonnegative matrix factorization //
Information Processing and Management: Volume 42 Issue 2, March 2006. – P. 373 – 386
3. Anatoly Anisimov, Oleksandr Marchenko, Andrey Nikonenko, Elena Porkhun, Volodymyr Taranukha: Ukrainian WordNet: Creation and Filling.
FQAS 2013. – P. 649–660
4. Tim Van de Cruys. Anon-negative tensor factor-ization model for selectional preference induction. Natural Language Engineering. – 2010. –
16(4). – P. 417–437.
5. Tim Van de Cruys, Laura Rimell, Thierry Poibeau, and Anna Korhonen. Multi-way Tensor Factorization for Unsupervised Lexical Acquisition.
In International Conference on Computational Lin-guistics (COLING), Mumbai, India, 08/12/2012-15/12/2012, P. 2703–2720. The COLING
2012 Organizing Committee. 2012.
6. Scott Deerwester, Susan T. Dumais, GeorgeW. Furnas, Thomas K. Landauer, and Richard Harshman. Indexing by latent semantic analysis // J.
OF THE AMERICAN SOCIETY FOR INFORMATION SCIENCE. – 1990. – 41(6). P. 391–407.
7. Brett W. Bader, Tamara G. Kolda, et al. Matlab tensor toolbox version 2.5. Available online, January. 2012.
8. Khushboo Kanjani. Parallel non negative matrix factorization for document clustering. May, 2007.
9. Volodymyr Kysenko, Karl Rupp, Oleksandr Marchenko, Siegfried Selberherr, and Anatoly Anisimov. Gpu-accelerated non-negative matrix
factorization for text mining. of Lecture Notes in Computer Science, Springer. – 2012. – V. 7337 – P. 158–163.
10. Chao Liu, Hung-chih Yang, Jinliang Fan, Li-Wei He, and Yi-Min Wang. Distributed nonnegative matrix factorization for web-scale dyadic data
anal-ysis on mapreduce. In Proceedings of the 19th In-ternational Conference on World Wide Web, WWW’10, New York, NY, USA. ACM. –
2010. – P. 681–690.
11. Daniel D. Lee and H. Sebastian Seung. Algorithms for non-negative matrix factorization. – 2000.
12. nVidia, 2013b. CUDA CUSPARSE Library. nVidia, August.
13. nVidia, 2013a. CUBLAS Library User Guide. nVidia, v 5.0 edition, October.
|