Оптимізація. Чорний ящик



У цій статті ми постараємося зрозуміти, яким чином пошук оптимальних параметрів об'єкта або моделі може бути автоматизований і зведений до оптимізації деякої ресурсномісткою функції типу «чорний ящик». Докладно розглянемо математичний метод (і його реалізацію) для ефективного вирішення цього завдання.

Трохи філософії

Якщо не вдаватися в деталі, то процес вирішення багатьох інженерних (і не тільки) завдань може бути зображений у вигляді такої схеми:



Простий приклад — проектування кузова спортивного автомобіля. В даному випадку дизайн — це геометрія кузова, модель — деякий розрахунковий пакет, а результат — скажімо, обчислюване аеродинамічний опір.

Важливо ось що — в наведеній вище схемі один із блоків завжди невідомий. І це приводить нас до трьох фундаментальних задач:

  • Пряма задача (невідомий результат)
  • Побудова моделі (невідома модель)
  • Обернена задача (невідомий дизайн)
Пряма завдання принципово проста. По суті це задача про побудову інтерфейсу — користувачеві повинно бути зручно задати дизайн та інтерпретувати результати роботи програми (подібні інтерфейси вже давно існують і продовжують удосконалюватися). У той же час побудова моделі — завдання принципово складна. Вчений, виходячи зі спостережень, розуміння і аналізу, займається науковою творчістю, яке в загальному випадку не автоматизируемо, принаймні на сьогоднішній день.

А ось із зворотним завданням все не так однозначно. З одного боку, процес створення дизайну також являє собою творчість. А з іншого боку, дизайн можна представити у вигляді набору деяких чисельних параметрів (найбільш істотних у розглянутому контексті), а результат — у вигляді деякої числової оцінки того, наскільки гарний дизайн. Таким чином шукається набір параметрів, які максимізують/мінімізують деяку величину, а значить зворотна задача — це задача математичної оптимізації і її рішення може і повинно бути автоматизовано.

до Речі кажучиПошук параметрів фіксованого моделі (наприклад, коефіцієнтів, що входять в її співвідношення) також можна розглядати як зворотну задачу. Тільки вирішується вона вже в контексті експерименту. У такому разі шукається набір параметрів, що мінімізують різницю між прогнозом моделі і результатом експерименту. Тобто завдання також зводиться до оптимізації.

І оскільки ми заговорили про оптимізацію, потрібно розуміти з якою функцією ми маємо справу. А функція ця непроста:

  • Ми знаємо тільки вхідні параметри (input) і, за результатами чисельного експерименту, значення скалярної оцінки (output). Сама ж функція невідома, оскільки являє собою складну систему. І тому простіше всього називати її чорним ящиком.
  • Ця функція може бути дуже ресурсномісткою і займати хвилини, години і дні для обчислення в одній точці (тобто для одного набору параметрів).
Таким чином, для розв'язання обернених задач нам потрібно навчитися ефективно оптимізувати такі чорні ящики (expensive black-box functions). Далі поговоримо про математику, яка для цього використовується.

Математичний принцип

Перед тим, як заглибитися в рішення, постараємося зрозуміти наступне. Головна проблема — це не те, що функція являє собою чорний ящик (не має аналітичного виразу), а те, що вона довго обчислюється. Якщо б вона обчислювалася швидко, то завдання її оптимізації вирішувалася яким-небудь простим проходом по сітці, методом Монте-Карло, або одним з градієнтних методів. Але, враховуючи дорожнечу функції, і відповідно дорожнечу обчислення її градієнта, потрібно діяти більш економно.

Зробимо так — спочатку обчислимо функцію в дуже невеликій кількості точок, рівномірно покривають простір пошуку, а потім интерполируем отримані значення деякої гладкою поверхнею. Таким чином ми реконструюємо вихідну функцію за значенням у кількох точках, не обчислюючи градієнти. Отримана поверхня вкаже нам на найбільш вірогідне місце розташування оптимуму, де ми вже зможемо його уточнити.

Початкова вибірка

Для простоти будемо вважати, що кожна змінна нашої функції має свій незалежний інтервал значень. Таким чином простір пошуку являє собою паралелепіпед в багатовимірному просторі. Якщо ж нормалізувати кожну змінну, то простір пошуку стає одиничним кубом (тут і далі будемо опускати префікс «гіпер»).

Перше питання — як найбільш ефективно (найбільш рівномірно) розташувати задане початкове кількість точок у просторі пошуку (кубі)?

В одновимірному випадку, коли функція задана на відрізку, рішення очевидно — рівномірна сітка. Однак рівномірна сітка працює помітно гірше в багатовимірному випадку оскільки кількість точок сильно прив'язане до розмірності простору. Скажімо, в чотиривимірному просторі сітка розміру 2 утворює image(16) точок, у той час як сітка розміру 3 утворює image(81) точку. При цьому, якщо у нас є обчислювальний час для, припустімо, 30 точок, 16 точок для нас — надто мало, а 81 — занадто багато. І ситуація, очевидно, буде збільшуватися із зростанням розмірності.

Ще один варіант — використовувати випадкову вибірку. У цьому випадку ми не прив'язані до розмірності простору і можемо генерувати довільну кількість точок. Однак, виникає інша проблема — нерівномірність покриття. Ось типове розташування 20 випадкових точок в одиничному квадраті:



На допомогу приходить метод під назвою Latin hypercube (LH). Уявімо, що потрібно розставити 8 не б'ють один одного тур на шахівниці — тобто в кожному рядку і кожному стовпці дошки має бути рівно одна човен. Це і є LH, тільки замість тур ми розміщуємо точки, а дошка може бути багатовимірної і її розмір може бути відмінним від 8. Можна показати, що кількість точок в цьому випадку дорівнює розміру дошки і не прив'язане до розмірності простору. При цьому існує безліч різних LH заданого розміру і нас будуть цікавити не все підряд, а тільки ті, в яких точки розташовані як можна більш рівномірно. Виявляється, що побудувати таке розташування можна за допомогою досить простого методу.

Для прикладу будемо розміщувати 20 точок у квадраті. Спочатку розташуємо точки діагональним чином:



Далі ми вибираємо 2 довільних рядка/стовпця і міняємо їх місцями. Після кожного такого обміну будемо стежити, щоб рівномірність покриття зросла. Зробити це можна, визначивши таку міру:

image

де image— евклідова відстань між двома точками, а image— кількість точок. Якщо після чергового обміну значення imageзменшилася (рівномірність зросла), ми зберігаємо це стан, в іншому випадку продовжуємо ітерації до тих пір, поки не досягнемо поліпшення. Якщо виконати такі ітерації деяке велике число разів (скажімо, 1000), то від діагонального розташування ми приходимо до наступної картині:



Це саме те, що нам потрібно — деякий заданий кількість точок, рівномірно заповнюють куб заданої розмірності.

Уточнюючі ітерації

Як тільки ми визначилися з початковою вибіркою, можна приступити до побудови интерполирующей поверхні для уточнюючих ітерацій. Будемо будувати поверхню за допомогою Radial basis functions (RBF) з кубічним базисом. Рівняння такої поверхні дуже просте і компактне:

image

Тут image— кубічний базис (image) imageimageimage— коефіцієнти, які знаходяться з умови інтерполяції та рішення відповідної лінійної системи рівнянь. Отримана поверхня imageпобудована з використанням imageточок imageі дозволяє передбачити значення вихідної функції в довільній точці image.

Після того, як початкова поверхня побудована, для уточнення оптимуму ми скористаємося ітераційним алгоритмом, детально описаним тут. Сенс алгоритму в тому, щоб знайти оптимум побудованої поверхні (можна навіть методом Монте-Карло, поверхня набагато «дешевше» вихідної функції) так, щоб нова точка знаходилась не ближче заданого мінімального відстані до всіх інших точок — щоб пошук не застрявав на одному місці. Це як якщо б ми розмістили в кожній точці обмежувач у вигляді м'ячика заданого радіуса, який би не дозволяв новій точці бути занадто близько. У методі вводиться функція щільності (відношення обсягу таких м'ячиків до об'єму куба):

image

де image— кількість уточнюючих ітерацій image— номер уточнюючої ітерації image— початкова щільність, а ступінь imageвизначає швидкість убування щільності з ітераціями.

Вибираючи ті або інші значення imageimage, можна задавати необхідний акцент на глобальному/локальному пошуку.

Приклад

Для демонстрації роботи методу розглянемо задачу мінімізації простої функції, наприклад, такий:

image

Будемо шукати її мінімум на області image, де вона виглядає ось так:



Удамо, ніби вираз функції нам невідомо (чорний ящик) і ми лише можемо обчислити її значення в зазначених точках. Запустимо наш метод, використовуючи 15 точок — 10 для LH і 5 для уточнюючих ітерацій:



Також запустимо метод, використовуючи 30 точок — 20 для LH і 10 для уточнюючих ітерацій:



Очевидно, що якість реконструкції вихідної функції з допомогою RBF-поверхні і точність методу безпосередньо залежать від заданого кількості точок.

Про глобальний оптимум і точності методу

Якби задана функція володіла деякими властивостями (по гладкості/опуклості/...), то напевно можна було б побудувати взаємозв'язок між кількістю точок і точністю методів. Однак наша функція — чорний ящик, а значить ми не знаємо про неї абсолютно нічого. Ніщо не заважає нашій функції виглядати ось так:



Цей штучний приклад наочно показує, що просто неможливо дати жодних гарантій щодо знаходження глобального оптимуму (а тим більше точності його знаходження).

Тому метод закладена дуже проста логіка — користувач замість заданої точності вказує доступний обчислювальний час (кількість точок), а метод всього лише ефективно використовує це обчислювальний час, намагаючись знайти глобальний оптимум.

Ядра, ядерця

Вище ми говорили про математичні методи, що дозволяють ефективно шукати оптимум, економлячи на обчисленнях функції. Але ми забули ще про одне — про паралельності. Чому б не обчислювати функцію в декількох місцях одночасно, використовуючи для кожного обчислення окреме ядро процесора? Адже ядер зараз чимало, а реалізувати таку фічу — справа нескладна. Ось так виглядає велосипед на Пітоні, який «мапит» задану функцію на кілька точок у паралельному режимі:

import multiprocessing as mp

def pmap(f,batch,n):

pool=mp.Pool(processes=n)
res=pool.map(f,batch)
pool.close()
pool.join()

return res

Ця проста штука дає приріст продуктивності, що дорівнює кількості використовуваних ядер!

Код

Код методу написана на Пітоні і швидше за все не є еталоном стилю, зате простий і працює. Зроблено все у вигляді одного модуля, що включає кілька функцій.

Посилання на GitHub

Від користувача вимагається надати чорний ящик у вигляді питоновской функції. Наприклад, у разі використання зовнішнього пакету (через системний виклик) вона може виглядати так:

def fun(par):

# modifying text file that contains design/model parameters
...
# performing simulation in external package
os.system(...)
# reading results, calculating output
...

return output

Аргументом функції має бути вектор параметрів, а результатом — невід'ємна скалярна величина. Запуск методу (для мінімізації функції) виглядає ось так:

from blackbox import *


def fun(par):
...
return output


if __name__ == '__main__':

search(

f=fun, # given function

resfile='output.csv', # .csv file where iterations will be saved

box=np.array([[-1.,1.],[-1.,1.]]), # range of values for each variable

cores=4, # number of cores to be used

n=8, # number of function calls on initial stage (global search)
it=8 # number of function calls on subsequent stage (local search)

)

Пара уточнень про те, як зберігаються і інтерпретуються результати, наведені в описі на GitHub.

На майбутнє

Хотілося б допилити дві речі:

  • Поки що аргументи функції задаються на незалежних інтервалах і простір пошуку являє собою куб. Якщо ж на аргументи накладені деякі обмеження, то простір пошуку буде мати довільну форму. Треба буде придумати, як найзручніше задавати такі обмеження в коді, і, що важливіше, як рівномірно розміщувати початкові точки в області довільної форми.
  • Непогано було б замінити RBF-поверхня на більш потужну (наприклад, імовірнісну). Це дало б деякий приріст до ефективності методу.


Література

Стаття з повним описом методу

Ще пара статей по темі:

Regis, R. G. Shoemaker and, C. A., 2005. Constrained global optimization of expensive black box functions using radial basis functions. Journal of Global optimization, 31(1), pp.153-171.

Holmström, K., 2008. An adaptive radial basis algorithm (ARBF) expensive for black-box global optimization. Journal of Global Optimization, 41(3), pp.447-464.

Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.