Численный анализ многокубитных систем

[С.И. Доронин]

Пипа, опять все правильно! Думаю, редукцию мы осилили .

[Pipa]

    Поясню причину первоначальных неурядец с редуцированием.
Как известно, позиционные системы счисления основаны на том, что число записывается в несколько разрядов, среди которых естественным образом выделяются младшие и старшие разряды.
    Если расположить натуральную последовательность чисел в ЛЮБОМ ПОРЯДКЕ (как в возрастающем, так и убывающем), то младший разряд легко узнаваем по его особенности меняться на каждом шаге, от числа к числу. Чем старше разряд, тем реже в нем происходят изменения. Например, в первой сотне чисел от 0 до 100, разряд единиц меняется каждый раз, разряд десятков меняется через десять шагов, разряд сотен сменится всего один раз.
    Теперь взглянем на то, как Доронин заполняет матрицу плотности:

                  <0|_A<0|_B   <0|_A<1|_B   <1|_A<0|_B   <1|_A<1|_B
                   ______________________________________
|0>_A|0>_B  |
|0>_A|1>_B  |
|1>_A|0>_B  |
|1>_A|1>_B  |
 

Для наглядности я распишу последовательность следования базисных векторов в таком виде:
A   B
0    0
0    1
1    0
1    1
Очевидно, что B выступает в роли МЛАДШЕГО разряда – ведь чередование в его ряду происходит чаще.
Еще более наглядно выглядит матрица для трехкубитной системы:
A   B   С
0    0    0
0    0    1
0    1    0
0    1    1
1    0    0
1    0    1
1    1    0
1    1    1
Как видим, добавление кубита C привело к тому, что младшим стал он.
Т.е. с ростом порядкового номера столбца и строки перебор базисных векторов строится таким образом, что чем ближе буква кубита к началу алфавита, тем старшее место он занимает.
     Используя битовую маску для редуцируемых кубитов, я по логике алгоритма должна была опираться на битово-цифровую систему, а не буквенно-алфавитную. И по началу считала кубиты в порядке возрастания старшенства.  Т.е. за кубит А я  принимала самый младший.
      Теперь, когда выяснилось, что порядок здесь обратный, я вывернула порядок наизнанку (произвела битриверс), после чего стали получаться правильные результаты.

[Pipa]

Доронин С.И.
    Устроила программе жесткую модерацию . Привела в порядок все, что можно, для того чтобы потом было меньше путаницы.
    Перво-наперво ввела классы матриц: RealMatrix, CompexMatrix, IntMatrix (для действительных, комплексных и целочисленных элементов, соответственно). И для полноты компании такие же векторы: RealVector, CompexVector, IntVector.
    Каждый из объектов этих классов аллокирует для себя память в конструкторе и освобождает ее в деструкторе. Хранят свои габариты. Содержат функцию Zer() для обнуления содержимого и другие полезные функции: нахождения следа, копировавания, распечатки в файл, чтения их файла и т.п.
    После этого программа стала совсем компактной и гораздо более понятной. Число параметров у функций уменьшилось, т.к. теперь можно передать в параметре ссылку на матрицу, не указывая ее размер. При комплексных матрицах получается один параметр вместо трех.
    Комплексная матрица CompexMatrix, включающая  в себя два массива (для действительной и мнимой части), образована на базе действительной RealMatrix. Это позволяет использовать ее в качестве параметра для функций, требующих действительную матрицу, и тогда будут использоваться только ее действительная часть.
    Операторы члены класса позволяют использовать индексацию элементов так же, как и было раньше. Подстановка возложена на компилятор, и исполнение программы это не задерживает.
   Например:
CompexMatrix ro(8,8), red(2,2);  // создание матриц 8x8 и 2x2
ro.Read( "in.txt");  // чтение матрицы из файла in.txt
Reduce( ro, red, "BC");  // редуцирование по кубитам B и C
red.Print( "out.txt");  // печать матрицы в файл out.txt
   Или
ro.re[2][3] = 5;  // присвоение значения действительной части элемента ro₂₃
ro.im[2][3] = -2;  // присвоение значения мнимой части элемента ro₂₃
   Или
x = ro.Trace();  // вычисление следа
   Или
ro.CopyTo( ro1);  // копирование матрицы ro в матрицу ro1
   и т.д.

    В ходе переработки программы возник вопрос, связанный с размерностями. Раньше этого не было заметно, но сейчас вышло на поверхность.  Класс BLOCK имеет две "габаритки". Словом "габаритки" я называют размер, определяющий длину массивов (векторов и матриц, членов класса BLOCK).
    1) первая габаритка N = 2^spin
    2) вторая габаритка cm = параметр конструктора, получающий значение из Hami->cm.
     При этом одни массивы имеют длину N, а другие cm, и от этого возникает формальная путаница. Например, в функции cosin() векторы размерностью N умножаются на матрицы размерностью cm x cm. Формально это нарушение.
     Еще хуже в функции BLOCK::Iz(). Там "крутятся" чужие массивы, принадлежащие HAMI (связь через параметры). Но циклы там по габариткам BLOCK'а, а не HAMI, который эти массивы аллокировал. Зачем это? Не лучше ли самому HAMI заниматься своими массивами, про размеры которых он в курсе?
     Размерности и в других местах имеют массу нестыковок, но программа продолжает работать только потому, что в ней N всегда совпадает с cm. Хоть и считает HAMI эту cm через lim, но точно так же, как и BLOCK вычисляет свое N.
     Это так и задумано было? Или все-таки могут быть режимы работы, когда cm < N?
     Нельзя ли чуточку пояснить относительно разбиения алгоритмов на BLOCK и HAMI? HAMI – это гамильтониан, с ним вроде бы все понятно. А что такое BLOCK? И может ли он быть больше гамильтониана?

P.S. Даже если в программе не всё потребуется для дальнейшей работы, все равно прошу пояснить эти моменты, т.к. хочется полностью распутать опасные перехлесты (грозящие выходом за пределов массивов), прежде чем развивать на ее базе дальнейшие вычисления.

[С.И. Доронин]

Пипа, спасибо! Исходники я тоже получил.

Размерности и в других местах имеют массу нестыковок, но программа продолжает работать только потому, что в ней N всегда совпадает с cm. Хоть и считает HAMI эту cm через lim, но точно так же, как и BLOCK вычисляет свое N.
     Это так и задумано было? Или все-таки могут быть режимы работы, когда cm < N?

Этот «хвост» тянется от задачи, которая решалась ранее. Для некоторых гамильтонианов, в частности для MQ-гамильтониана, который сейчас задействован, можно перейти к задаче меньшей размерности – в одном из базисов матрица гамильтониана и сама матрица плотности разбивается на блоки, и для полного решения задачи достаточно работать с одним блоком, матрица которого в два раза меньше исходной. Поэтому и был предусмотрен режим, когда cm < N.
Но нам, я думаю, лучше работать с полными матрицами, поскольку для других гамильтонианов этого не будет. Поэтому, cm лучше везде заменить на N, чтобы не путаться с размерностями.

Нельзя ли чуточку пояснить относительно разбиения алгоритмов на BLOCK и HAMI? HAMI – это гамильтониан, с ним вроде бы все понятно. А что такое BLOCK? И может ли он быть больше гамильтониана?

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

[С.И. Доронин]

Пипа, я посмотрел Ваш последний релиз программы. Возник такой вопрос – трехкубитная система у Вас считается? У меня программа вылетает на функции Tqli с сообщением "Too many iterations in tqli". Причем это только для 3 кубитов, с другим числом кубитов считает нормально. В предыдущем варианте программы, который Вы высылали перед этим, такой проблемы не было, трехкубитная система считалась без проблем.

[Pipa]

   Трехкубитную систему ровно с тем же сообщенем "Too many iterations in tqli" не считала исходная ваша программа, откомпилированная мной и минимальными изменениями только (вставка дополнительных хидеров) на компиляторе Boraland C++ 5.02.
   Отчего мой прошлый вариант посчитал 3-кубита это удивительно и для меня. По крайней мере, сравнить результат мне не с чем, поскольку вместе и исходниками, работающего exe-модуля вы не прислали. (А нельзя ли прислать в качестве "эталона"?)
   С причиной ошибки разберусь, а интуитивно мне кажется виновата плохая проверка малости fabs(e[m]+dd)==dd
   Например, Уилкинсон в своем знаменитом Справочнике, рекомендует в этом месте условие: fabs(e[m]) <= DBL_EPSILON*dd

[Pipa]

С.И. Доронин
    Мое первоначальное предположение оказалось неверным – в вычислительном плане оба условия на расщепление тридиагональной матрицы являются эквивалентными. Но то, что сбой обнаружился в процедуре QL-алгоритма меня сильно ... порадовало .
    Порадовало оттого, что моя любовь к матрицам воспылала от прочтения замечательной книги Парлетта (B.N. Parlett). Какая не скажу , тем более что в русском переводе у него вышла всего одна книжка, которая нынче представляет библиографическую редкость. Вообще-то практической пользы от той книги вроде бы и нет, т.к. проблема собственных значений для симметричных матриц была решена давно (около 1960 г) и с тех пор ею мало кто интересуется, т.к. все эти алгоритмы давно стали достоянием любого библиотечного пакета для матричных операций (EISPACK, LAPACK, SCILAB и многих других). Но, тем не менее, книга написана настолько хорошо, что ее чтение доставляет эстетическое удовольствие. В ней даже в помине нет ничего педагогически-занудного. И ничего популярного тоже нет! Вот ведь как бывает.
     Книжку вспомнила оттого, что в ней были затронуты вопросы, относящиеся к сходимости QL-алгоритма. Когда устанавливали порог в 30 итераций на одно собственное значение было даже заключено пари на кругленькую сумму , что никто не найдет такой матрицы, которая достигнет этого лимита (обычно более 4-х бывает редко). И хотя сроки того пари давно истекли (нет уже в живых того, что бился об заклад), но я была крайне возбуждена, увидев, что 3-кубитную систему вышибло при переходе за тридцатник.
     К сожалению, прославить свое имя такой удивительной матрицей мне не пришлось , причина несходимости итераций оказалась гораздо более прозаической – машинные нули – вещь давно известная.
     И, тем не менее, матрица гамильтониана меня удивила – вопреки моим ожиданиям она имела практически нулевой след (e-307), а крупные числа (e+4) оказались вне диагонали. Из-за симметрии матрицы гамильтониана, его приведение к тридиагональной форме сохраняет нули на главной диагонали. Да ладно бы нули, а там шушера с порядками от e-309 до e-290, т.е. на самом краю процессорной арифметики. К тому же, неудачная процедура масштабирования в tred2() еще более усугубляла ситуацию.
     Процедуру «тридиагонализации» (tred2) пришлось видоизменить – привести ее к виду одноименной процедуры из пакета EISPACK. Сбой от этого пропал, но сильно радоваться этому не стоит – условия сходимости одинаково плохи для данного типа гамильтониана при ЛЮБОМ числе кубитов. Остальным просто повезло, что они уложились в 30. А в 3-кубитной системе балансировало на волоске, из-за этого то сходилось, а то нет.
     Ранее мне приходилось встречаться с задачей нахождения собственных значений, но мои матрицы всякий раз были типа корелляционных/ковариационных, т.е. явно тяготели к диагонали (именно на ней были или единицы или самые крупные числа). Соответственно этому, таких проблем с ними никогда не было.
     Знаю, как составляется гамильтониан для N-спиновой системы в ПМР. Там гамильтониан блочный:
2² = 4 = 1+2+1
2³ = 8 = 1+3+3+1
2⁴ = 16 = 1+4+6+4+1
2⁵ = 32 = 1+5+10+10+5+1
и тоже с тяготением к диагонали (даже углы у него голые). А этот гамильтониан такой из себя странный...
     В связи с этим хотелось бы попросить подробнее рассказать про данный тип гамильтониана. Тем более что вы ранее уже обмолвились, что у вас имеется несколько способов их построения.

[С.И. Доронин]

В конце мая в Письмах в ЖЭТФ (самый «крутой» физический журнал у нас в России ) вышла наша лабораторная статья: http://www.jetpletters.ac.ru/ps/1346/article_20329.shtml
С.И. Доронин, А.Н. Пырков, Э.Б. Фельдман, Запутанность в альтернированных открытых цепочках ядерных спинов s=1/2 с XY-гамильтонианом, Письма в ЖЭТФ 85, 627 (2007).

Она на русском языке и, по-моему, свободна для скачивания. Это не совсем то, о чем мы говорим здесь, но все равно близко. Там редуцированная матрица плотности (только двухкубитная) найдена аналитически для любого числа спинов в цепочке. Правда и взаимодействия кубитов учитывались не все, а лишь ближайших соседей (XY-гамильтониан).

P.S. Пипа, я извиняюсь, что не ответил на Ваш вопрос насчет гамильтониана, разговор тут долгий. Кратко могу сказать, что в той задаче, которая сейчас у нас «крутится» используется так называемый многоквантовый (МК) гамильтониан. В отличие от обычного гамильтониана в ЯМР, который отвечает за квантовые переходы лишь между соседними энергетическими уровнями, МК-гамильтониан реализует все возможные переходы между квантовыми уровнями, он используется в многоквантовом ЯМР.

[Pipa]

Oleg Shaposhnikov
   Прошу меня извинить, но заведенный вами разговор мне пришлось выделить в отдельную тему "Гадание на кубитах" http://quantmag.ppole.ru/index.php?option=com_smf&Itemid=34&topic=180.0
   Вы вправе изменить ее название, отредактировав свое первое сообщение в ней.

Модератор раздела

[ksv]

Сергей Иванович, Пипа скажите, а в каком состоянии сейчас находится сей проект?
Что уже сделано, что еще планируется сделать?

[С.И. Доронин]

ksv

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


Это динамика квантовой запутанности на различных квантовых уровнях системы из 4 кубитов, связанных диполь-дипольным взаимодействием (здесь используется обычный дипольный гамильтониан, а не многоквантовый) при начальном условии |Ψ>=|0001>. Т.е. в начальный момент времени система находилась в сепарабельном состоянии, но в результате взаимодействия возникает запутанность, и мы можем видеть ее динамику.

Например, запутанность одного кубита А с остальными кубитами, т.е. запутанность типа A–BCD имеет вид:


Запутанность двух кубитов AB с другими два (CD), эволюционирует следующим образом:


Это второй квантовый уровень.

На третьем, самом «тонком» квантовом уровне четырехкубитной системы, т.е. запутанность типа ABC–D ведет себя так:


Самая простая двухкубитная запутанность типа A–B имеет такую динамику:


Хочу поблагодарить Пипу, – она просто молодец! Это ее усилиями программа появилась на свет, и основной движок, на данный момент, практически полностью готов и протестирован. Сейчас продолжается работа над графическим интерфейсом, чтобы более удобно ей было пользоваться. В дальнейших планах – наращивание числа различных используемых гамильтонианов, ну и различная мелочевка…

Хорошо, если Пипа поделится своими впечатлениями от работы над программой, думаю, ей есть что сказать , хотя бы общую оценку…

Я полагаю, делается очень большое и нужное дело, уже сейчас ясно, что мы получаем в свои руки чрезвычайно мощный и очень эффективный инструмент по анализу квантовой запутанности, причем в ее динамике. Программа многое может дать в плане понимания сути квантовой запутанности и механизма ее действия. Причем, уже на стадии отладки и тестирования, получен ряд интересных и важных теоретических выводов, которые я потихоньку осмысливаю . То ли еще будет…

[ksv]

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

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

[Pipa]

Насколько понимаю динамика запутанности сильно зависит от гамильтониана.

    Не то слово! Но было еще более странно, если бы от гамильтониана ничего не зависело. В конце концов, гамильтониан и задает динамику.

А гамильтониан задаете матрицей? Как она выглядит?

     Сергей Иванович, кажется, посылал вам исходник той програмы, от которой начали танцевать. Там есть и алгоритм постороения гамильтониана. Лишь совсем недавно Доронин предложил включить и другой гамильтониан. Причем оба они задаются не готовой матрицей, а некоторым алгоритмом (правилом), который строит такой гамильтониан для заданого числа кубитов. Оба гамильтониана представляют собой действительные симметричные матрицы, причем достаточно разряженные.
     Для малого числа кубитов они выглядят примитивно, например для двух кубитов первый выглядит так:
0       0       0       1
0       0       0       0
0       0       0       0
1       0       0       0
а второй так:
1       0       0       0
0      -1      -1       0
0      -1      -1       0
0       0       0       1
(нормировочный множитель опущен).
    И все-таки очень странно, что запутанность очень сильно зависит от матрицы исходного состояния. Мне раньше казалось, что процесс должен "окончательно запутаться", т.е. придти к некоторому равновесному состоянию. Однако этого не происходит. Эволюция сводится к процессу "вращения состояния" в многомерном пространстве, причем положение "оси вращения" определяет структура гамильтониана. Вращаться эта "штуковина" может долго, и видимо можно дождаться, когда система вновь пройдет исходное состояние. Т.е. процесс носит явно периодический характер, несмотря на то, что период может быть очень большим (ожидание, когда множество осцилляторов вновь окажется в одинаковой фазе).
     Открытое обсуждение этой темы с некоторых пор было прервано по просьбе Доронина, из-за того, что предложенный им метод расчета запутанности до сих пор еще не опубликован в открытой печати. С сегодняшнего дня он находится в отпуске. А моя интенсивная работа в последнее время была обусловлена необходимостью сотворить на скорую руку такой програмный инструмент, чтобы во время отпуска у него была возможность более тщательно исследовать закономерности динамики запутанности (ноутбук у его с собой будет, а интернет нет) . Вследствие этого, усилия были направлены по большей мере на интерфейс пользователя, чтобы дать возможность задавать условия задачи без программирования, а исключительно средствами пользовательского интерфейса.
    В настоящий момент программа представляет не готовый продукт, расчитанный на широкий круг пользователей, а скорее макет, в котором все происходящие события по максимому "оголены". Т.е. сейчас это в высшей степени отладочный инструмент, с максимальной прозрачностью и средствами просмотра промежуточных результатов. К этому вынуждает сама вычислительная сложность задачи, не позволяющая "на глазок" оценивать правдопобность результатов, что заставляет выводить многие дополнительные характеристики промежуточных продуктов в виде текстовых таблиц.
    Что же касается программирования в новой для меня системе Borland C++ Builder 6.0, то я была просто удручена тем, насколько это занятие не вяжется с тем, что я считала программированием раньше. Тут играет основную роль не знание языка программирования, а знание того, на какие кнопочки надо нажимать, и где их искать. Такое занятие больше похоже на игру Sims , чем на программирование, - выбираешь мышкой из списка тело героя, голову, цвет волос, прическу, одежду и т.д., вплоть до меблировки комнаты. Не программирование, а какая-то игра типа "Одень Барби" . Вот и маешься, перетрясая бесконечный гардероб в поисках того, что сгодилось бы для программы.

[Pipa]

В конце мая в Письмах в ЖЭТФ (самый «крутой» физический журнал у нас в России ) вышла наша лабораторная статья: http://www.jetpletters.ac.ru/ps/1346/article_20329.shtml
С.И. Доронин, А.Н. Пырков, Э.Б. Фельдман, Запутанность в альтернированных открытых цепочках ядерных спинов s=1/2 с XY-гамильтонианом, Письма в ЖЭТФ 85, 627 (2007).
Она на русском языке и, по-моему, свободна для скачивания. Это не совсем то, о чем мы говорим здесь, но все равно близко.

    Если что и действительно близко к тому о чем здесь говорят, то это, на мой взгляд, статья Городецких:
Horodecki M, Horodecki P. Reduction criterion of separability and limits for a class of protocols of entanglement distillation,
увидевшая свет почти 10 лет тому назад.
    На эту статью набрела совершенно случайно – ее мне выдал поисковик Google, когда я искала совершенно иной материал. Однако слово «entanglement» (запутанность) привлекло мое внимание и заставило взглянуть на текст этой статьи, благо это был препринт со свободным доступом для скачивания.
    Проглядев статью «по диагонали», у меня появился страшный зуд обвинить Доронина в плагиате . Настолько близким показался мне этот материал к тому, что сейчас делается.
    В своих работах Доронин не раз ссылался на «клон» Городецких, которые, по-видимому, в этой области являются заслуженными светилами. Однако все те ссылки относятся к «трехголовому клону», где присутствует еще Horodecki R, которого в этой статье нет.
    Не имея прямого интернет-доступа к Phys. Rev. Lett., я не cмогла прочесть более поздние статьи чрезвычайно плодовитого «клона» Городецких. Очень возможно, что в других своих статьях они повторяют, а то и развивают те идеи, которые я обнаружила в этом препринте 1998 года. И, тем не менее, даю ссылку на эту работу. Она, несомненно, того заслужила.

[С.И. Доронин]

Pipa

Естественно, что мера запутанности, которая используется в программе, возникла не на пустом месте . Помимо той статьи, которую Вы указываете M. Horodecki and P. Horodecki, Phys. Rev. A 59, 4206 (1999), можно назвать еще один первоисточник:
N.J. Cerf, C. Adami, and R.M. Gingrich, Phys. Rev. A 60, 898 (1999), в свободном доступе http://arxiv.org/abs/quant-ph/9710001
На эти статьи я ссылаюсь, только нужно понимать разницу между критерием сепарабельности (в данном случае речь идет о критерии приведения - reduction criterion) и мерой запутанности. Критерии сепарабельности обычно формулируются в виде неравенств и дают лишь общий ответ, является состояние сепарабельным или нет. Они не дают количественной величины (меры) запутанности. В свою очередь, на основе критериев сепарабельности, можно вводить и меры запутанности.

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

Чтобы лучше понять ситуацию могу привести такой пример. На основе РРТ-критерия (в файловом архиве есть статьи по этому критерию) мера запутанности была введена только спустя несколько лет в работе G. Vidal and R.F. Werner, Phys. Rev. A 65, 032314 (2002) http://arxiv.org/abs/quant-ph/0102117 . И это отдельная статья, в которой вводится и обосновывается именно мера запутанности (negativity) на основе ранее хорошо известного РРТ-критерия. Но эта мера не очень удобная в практическом применении, там есть ряд вопросов, которые затрудняют ее применение в той задаче, которая решается у нас.