[Oldring 0]

Узкое место в скорости расчетов - вычисление расширенной матрицы. Алгоритм простой, но сложность почти кубическая.

 

Там N*(N+1)/2 скалярных произведений векторов. Скалярные произведения обычно считаются на современных процессорах эффективно, без излишних накладных расходов, с использованием MAC или вообще SIMD команд. Но может потребоваться кодирование этой операции на ассемблере.

 

Прошу, приведите пример. Не понял, что подразумевается под устойчивостью.

 

http://ru.wikipedia.org/wiki/%D0%92%D1%8B%...%81%D1%82%D1%8C

[getch 0]

Для такой задачи естественно нормировать вектор коэффициентов так, чтобы Σx_i=1

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

При такой постановке задачи (ноль) количество решений бесконечно (это на первый взгляд).

 

Видится логичным такая ситуация:

Допустим, мы нашли искомый вектор (Vector), в котором дисперсия вектора (Matrix * Vector) минимальна.

Если мы меняем знак i-го столбца в Matrix, то новое решение было бы тем же вектором Vector, но только у которого i-й элемент тоже поменял знак.

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

Изменено 8 сентября, 2010 пользователем getch

[Oldring 0]

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

 

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

 

Правильное ограничение обычно получается из физики задачи. Расскажите, зачем вам это нужно - можно будет что-то подсказать.

[getch 0]

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

Конечно, будет ноль. Я только не понял, почему условие равенства единице первого элемента искомого вектора - это нехорошо. Вроде, при таком условии и обозначенное выше свойство решения будет выполняться.

Правильное ограничение обычно получается из физики задачи. Расскажите, зачем вам это нужно - можно будет что-то подсказать.

Есть такая наука - эконометрика. Стал с ней знакомиться и столкнулся с моей точки зрения с псевдо-научными рассуждениями и результатами (там не по делу (мое мнение) применяются довольно продвинутые результаты теорвера и матстатистики) эконометрики, за которые иногда авторы получали даже Нобелевские премии. Хочу разобраться просто, что откуда ведет. Благо огромное количество истории различных экономических показателей имеется. Разрабатываю инструментарий для исследования огромных баз экономических данных. Почему экономика, а не нечто другое? Причина только одна, нигде не накоплено такое большое количество информации. Ну нет измеренных данных изменения численности муравьев в муравейнике (и его роста) и подобных. В экономике есть.

Зачем мне это надо? Чтобы понимать или иметь хотя бы минимальное представление.

Изменено 8 сентября, 2010 пользователем getch

[Oldring 0]

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

 

Потому что минимум будет другим. Выделение первого вектора при этом нарушает симметрию между веаторами.

 

Если вам нужно обязательно чтобы при изменении знака вектора просто менялся знак его коэффициента - возьмите в качестве ограничения равенство суммы квадратов коэффициентов единице. Оптимальный вектор коэффициентов при этом будет собственным вектором, соответствующим минимальному собственному значению корреляционной матрицы.

 

Про эконометрику не скажу ничего, какое там ограничение "естественно". Скорее всего, там много вымысла и приближенности. Сложно ожидать линейности от экономики. :)

[getch 0]

Потому что минимум будет другим. Выделение первого вектора при этом нарушает симметрию между веаторами.

 

Если вам нужно обязательно чтобы при изменении знака вектора просто менялся знак его коэффициента - возьмите в качестве ограничения равенство суммы квадратов коэффициентов единице. Оптимальный вектор коэффициентов при этом будет собственным вектором, соответствующим минимальному собственному значению корреляционной матрицы.

Допустим, первый коэффициент равен Const1, и мы находим решение Vector1.

Если первый коэффициент равен Const2, находим решение Vector2.

При этом есть замечательное свойство: Vector2 = Vector1 * Const2 / Const1.

 

Поэтому решаю задачу так:

Первый коэффициент делаю единицей и нахожу решение Vector. Затем домножаю Vector на коэффициент, чтобы сумма квадратов его элементов была единица. Тем самым нахожу однозначное симметричное решение симметричной задачи.

 

Или я где-то что-то сильно упускаю?

Изменено 8 сентября, 2010 пользователем getch

[Oldring 0]

Или я где-то что-то сильно упускаю?

 

Да. Сделайте константой второй коэффициеинт. Получите в общем случае другой оптимальный вектор. Непропорциональный первому.

[getch 0]

Да. Сделайте константой второй коэффициеинт. Получите в общем случае другой оптимальный вектор. Непропорциональный первому.

Действительно, получаются непропорциональные решения.

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

И система уравнений получается вовсе не линейная. Тут, видимо, надо использовать какой-то численный метод решения. МНК подходит для этого?

 

Да. Сделайте константой второй коэффициеинт. Получите в общем случае другой оптимальный вектор. Непропорциональный первому.

И все же совершенно не понимаю природу этого.

Делаю так:

Определеляю первый коэффициент единицей, нахожу оптимальный вектор Vector1.

Заново решаю задачу, но уже определяю константой второй коэффициент, равным второму элементу из Vector1. Получаю оптимальный вектор Vector2.

Почему Vector1 не равен Vector2?!

Похоже, понял причину.

Изменено 9 сентября, 2010 пользователем getch

[getch 0]

Вроде, это задача квадратичного программирования. Пока смотрю материал на тему седловой точки функции Лагранжа. Загруз еще тот... Определиться бы с направлением.

 

Напоминаю формулировку задачи:

Найти такой вектор V, чтобы дисперсия вектора (Matrix*V) была минимальна.

При этом сумма КВАДРАТОВ элементов вектора V равна единице.

Мат. ожидание столбцов матрицы Matrix равно нулю.

[Oldring 0]

Действительно, получаются непропорциональные решения.

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

И система уравнений получается вовсе не линейная. Тут, видимо, надо использовать какой-то численный метод решения. МНК подходит для этого?

 

Нет, не подходит. С квадратичным ограничением задача сразу становится иная. Но также полностью обсосанная. Напишу еще раз, что решение есть собственный вектор, соответствующий минимальному собственному значению корреляционной матрицы. Собственные значения с собственными векторами - это также очень важные понятия из линейной алгебры, для их поиска разработаны свои эффективные и устойчивые алгоритмы, которые по своей природе итеративные, в отличие от решения системы линейных уравнений. В Матлабе их находит функция eig, для матрицы 10х10 она работает мгновенно. Но можете искать минимум квадратичной формы при квадратичном ограничении и явно. Или, например, вот таким способом, работающим для матриц с не слишком большыми числами обусловленностями:

 

1. Вычислить корреляционную матрицу сигналов: A=V'*V;

2. Найти обратную матрицу: B=inv(A);

3. Несколько раз возвести обратную матрицу в квадрат, не забывая её при этом перенормировать: B=B/max(max(abs(B )));B=B*B;

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

4. Найти столбец x матрицы B с наибольшей суммой квадратов элементов. Перенормировать его на равенство суммы квадратов единице: x=x/sqrt(x'*x). Найденный x и есть искомый столбец коэффициентов, минимизирующий дисперсию.

[getch 0]

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

Корреляционная матрица составляется, как корреляция между соответствующими столбцами входной матрицы?

 

Но можете искать минимум квадратичной формы при квадратичном ограничении и явно.

Под явный поиск что имеете в виду?

 

Или, например, вот таким способом, работающим для матриц с не слишком большыми числами обусловленностями:

 

1. Вычислить корреляционную матрицу сигналов: A=V'*V;

2. Найти обратную матрицу: B=inv(A);

3. Несколько раз возвести обратную матрицу в квадрат, не забывая её при этом перенормировать: B=B/max(max(abs(B )));B=B*B;

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

4. Найти столбец x матрицы B с наибольшей суммой квадратов элементов. Перенормировать его на равенство суммы квадратов единице: x=x/sqrt(x'*x). Найденный x и есть искомый столбец коэффициентов, минимизирующий дисперсию.

Что здесь обозначает V и abs? Перенормировка B на каждой итерации?

[Oldring 0]

Корреляционная матрица составляется, как корреляция между соответствующими столбцами входной матрицы?

 

Да. Произведение V'*V

 

Что здесь обозначает V и abs? Перенормировка B на каждой итерации?

 

Да. abs - абсолютное значение элементов матрицы. Точный вид нормировки не важен.

[getch 0]

3. Несколько раз возвести обратную матрицу в квадрат, не забывая её при этом перенормировать: B=B/max(max(abs(B )));B=B*B;

Дважды max - это опечатка?

[Oldring 0]

Дважды max - это опечатка?

 

Это функции Матлаба.

[getch 0]

1. Вычислить корреляционную матрицу сигналов: A=V'*V;

2. Найти обратную матрицу: B=inv(A);

3. Несколько раз возвести обратную матрицу в квадрат, не забывая её при этом перенормировать: B=B/max(max(abs(B )));B=B*B;

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

4. Найти столбец x матрицы B с наибольшей суммой квадратов элементов. Перенормировать его на равенство суммы квадратов единице: x=x/sqrt(x'*x). Найденный x и есть искомый столбец коэффициентов, минимизирующий дисперсию.

Воспроизвел, спасибо!. Ресурсоемкость минимальная при таком подходе.

Не знаю, правда, является ли найденный вектор оптимальным.

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

Подскажите, где можно почитать на эту тему? Или это все тот же MIT-курс лекций освещает?

 

3. Несколько раз возвести обратную матрицу в квадрат, не забывая её при этом перенормировать: B=B/max(max(abs(B )));B=B*B;

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

Мне хватает и 5-и итераций (корреляционная матрица 6x6). При такой быстрой сходимости норму разности между итерациями можно использовать простейшую: сумма квадратов элементов.

Изменено 9 сентября, 2010 пользователем getch