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

[Pipa]

С.И. Доронину

    В тексте программы есть функция MultiIx(), которая никогда не вызывается. Судя по содержанию, она считает экспоненты, причем не комплексные, а исключительно действительные.
    Нельзя прояснить вопрос о ее назначении?
   

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

Pipa

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

[Pipa]

С.И. Доронину

    Наконец-то дошли руки откомпилировать и запустить вашу программу. С компиляцией проблем не возникло (компилятор: Borland C++ 5.02).
    Основное время действительно пожирает функция Multi(), занимающаяся перемножением комплексных матриц. На этом шаге скорость можно увеличить, если отказаться от вычисления скалярных произведений каждой пары векторов в комплексной форме, а перейти на действительную. Сделать это можно так:
    Если r1[] и i1[] – соответственно действительная и мнимая части первого вектора, а r2[] и i2[] - второго, то их комплексное скалярное произведение может быть вычислено как сумма:
    sum_re = Summa(r1[]*r2[]) - Summa(i1[]*i2[])
    sum_im = Summa(r1[]*i2[]) + Summa(r2[]*i1[])
где все четыре суммы вычисляются по соответственным парам обычных double-векторов, что можно вычислить значительно быстрее (способ я подробно описала ранее в #7), не связываясь с комплексным представлением.
    К сожалению, для этого надо транспонировать матрицы так, чтобы вектора умножались не по столбцам, а по строкам. А это существенная переделка, в ходе которой легко допустить ошибку. Необходимость состоит в том, что в матрицах данного типа СТРОКИ содержат элементы, расположенные физически друг за другом, и переход на следующую пару элементов осуществляется простым инкрементированием адреса, когда как перебор по столбцу обходится значительно дороже.
    Нынешний вариант допускает напрасную трату времени еще и на формирование комплексной формы из двух ее частей на каждом шаге цикла. Это непроизводительно.

[Pipa]

С.И. Доронину

Есть 9-кубитная система!

    Оптимизация оказалась сильно проще, чем я думала – транспонировать ничего не пришлось, поскольку матрица cosin оказалась симметричной, а, значит, могла быть использована как уже транспонированная.
    Отказ от комплексной арифметики и ассемблерное вычисление векторных скалярных произведений дало ускорение в 3 раза (точнее в 3.2 - 3.3). Выходные файлы ускоренной и исходной программы полностью совпали для 5,6,7,8-кубитных систем (9-кубитная на исходной программе не запускалась из-за долгого ожидания результатов).
     В программу добавлена индикация времени исполнения. Получены следующие значения:

Кубитов        Исходная        Оптимизированная
      5                    1 сек               1 сек
      6                    8 сек               3 сек
      7                  56 сек             17 сек
      8                448 сек           134 сек
      9                    ?                1118 сек ( 18.6 мин)

Итого, менее 20 мин для 9-кубитной системы на компьютере Pentium-4 3.2 GHz с шиной 400 MHz.

    В вычислительном плане рабочие массивы rab1 и rab2 не нужны, поскольку результат работы функции Multi() утилизируется в матрице rm[][]. Вычисление ее элементов возможно сразу в общем цикле (вынесение ее вычисления в отдельный цикл избыточно, из-за чего и возникли рабочие массивы).
     В методическом плане меня настораживает, что элементы rm[][] вычисляются как сумма действительной и мнимой частей комплексного числа. "По духу" на это место просится модуль комплексного числа, т.е. корень из суммы квадратов реальной и мнимой части. Вы уверены, что там должна быть сумма, а не модуль?

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

Отлично! Пипа, Вы молодец!  Около двух минут для 8-кубитов это терпимо. На втором этапе, когда будет считаться мера квантовой запутанности, время тоже много уйдет. Там собственные значения надо будет вычислять на каждом шаге по времени, причем из комплексной эрмитовой матрицы, насчет метода тоже здесь надо думать, правда собственные вектора уже не нужны будут, что облегчает задачу.

В методическом плане меня настораживает, что элементы rm[][] вычисляются как сумма действительной и мнимой частей комплексного числа. "По духу" на это место просится модуль комплексного числа, т.е. корень из суммы квадратов реальной и мнимой части. Вы уверены, что там должна быть сумма, а не модуль?

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

Нам нужно будет вещественную и мнимую часть матрицы плотности хранить в отдельных файлах.
Да, я еще упоминал, что вместо вектора, который засылается в функцию Multi() в качестве параметра у нас должна быть матрица, т.е. умножение немного усложнится. Пока можно тот вектор, который был раньше, расположить на диагонали новой матрицы.

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

ksv

Но я говорил о другом. Сама форма уравнения, которае задает динамику системы. По моему мнению уравнение Лиувилля в этоим плане принципиально отличается от приводившейся мной системы уравнени, когда результат вычисления первого уравнения используется для расчета во втором уравнении, результат которого в свою очередь используется в первом уравнении на следующем шаге и т.д. В последнем случае время "возникает" в результате взаимодействия подсистем и их совместной эволюции, а не задается изначально в виде параметра. Другими словами, в уравнении Лиувилля мы можем рассчитать матрицу плотности системы на произвольном шаге t=a (где а - сколь угодно далеко отстоит от t=0), минуя все промежуточные расчеты матриц плотности на шагах от t=1 до t=a-1 путем применения оператора сдвига во времени. А в моем примере этого сделать принципиально не получится: нельзя рассчитать характеристики системы на шаге t=a без расчета всех предыдущих циклов взаимодействия в моменты 0<t<a, от которых зависит состояние системы в момент t.

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

Вот только не понимаю как в матрице плотности можно выделить отдельный кубит? Матрица ведь представляет общее (интерференционное) состояние всех кубитов системы.

Матрицы плотности для отдельных кубитов и более сложных подсистем типа AB, ABC, и т.д. находятся взятием частичного следа по окружению, получаются редуцированные матрицы плотности уже для подсистем.

А что вообще дает полный анализ квантовой запутанности? Т.е. насколько понимаю мера квантовой запутанности - это вещественный коэффициент от 0 до 1, показывающий степень взаимосвязи частей системы: от 0 (полностью сепарабельная система) до 1 (чистое состояние). Какие выводы можно сделать по динамике этих мер запутанности? Почему им придают такое большое значение?

Да, смысл запутанности примерно такой, но не совсем, максимальная мера запутанности это не просто чистое состояние, поскольку и сепарабельное состояние тоже может быть чистым, напр. |ψ>=|000…0>. Мера запутанности характеризует нелокальность, ее значимость, роль и «вес» квантовых корреляций по сравнению с классическими взаимодействиями. Мера запутанности, равная 1, означает, что в этом случае нет классических взаимодействий и локальных объектов. Мы можем посчитать несколько таких «тонкоматериальных» уровней, для которых существенны нелокальные корреляции, и посмотреть, как они взаимосвязаны друг с другом.

Как я понимаю смотреть надо только на алгоритм в функции Multi() ?

В этой функции только умножение матриц и получается нужная нам матрица плотности. Особо копать программу нет необходимости, достаточно иметь хотя бы общее представление. Как говориться, свою программу порой написать легче, чем разобраться в чужой . А Вы сможете написать какую-нибудь функцию? Например, взятие частичного следа и получение редуцированной матрицы плотности. Общий алгоритм есть в статье, о которой я уже упоминал http://quantmagic.narod.ru/volumes/VOL112004/abs1123.html на 5 странице.

[Pipa]

Мне нужна была сумма, там считаются интенсивности многоквантовых когерентностей, и правильность расчетов многократно проверена другими методами.

    А если так, то вы делаете двойную работу. Объясняю.
Раз уж вы эти две суммы все равно потом складываете

    sum_re = Summa(m2[k]*cos[k] - rm1[k]*sin[k])
    sum_im = Summa(m2[k]*sin[k] + rm1[k]*cos[k])
    rm[j] = sum_re + sum_im

то зачем же тогда их считать отдельно?
Накапливайте тогда их в одном сумматоре:
   
    rm[j] = Summa(m2[k]*cos[k] - rm1[k]*sin[k] + m2[k]*sin[k] + rm1[k]*cos[k])

и тогда увидите, что одно и тоже число умножается на cos и на sin, а потом эти произведения складываются. Отсюда следует, что если использовать вместо таблиц cos[] и sin[], таблицы их сумм и разностей cps[] = cos[]+sin[] и cms[] = cos[]-sin[], то матрица плотности может быть вычислена вдвое быстрее:

    rm[j] = Summa(m2[k]*cms[k] + rm1[k]*cps[k])

Т.е. если задача состоит в нахождении суммы действительной и мнимой частей сумм, то просто глупо "честно" множить матрицы. А поскольку эта процедура лимитирует скорость всей задачи, то такие действия уже просто расточительны.

[Pipa]

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

    В чем же дело? Пусть тогда другие задачи не пользуются этой фукцией Тем более когда универсальность достигается такой дорогой ценой.
    О чем базар? Перемножение матриц - это всего три строки кода! Из которых первые две - вложенные друг в друга циклы для обхода всех клеток результата, а третья вычисляет скалярное произведение пары векторов. Комплексный случай сложнее только тем, что третья строка у него чуточку длиннее (там уже не одно скалярное произведение, а сумма или разность четырех).
    Что мешает иметь функцию для перемножения матриц, когда нужны комплексные элементы произведения, и отдельно функцию, когда мнимая часть рассматривается, как действительная, и оттого может быть с ней сложена в одно действительное число? Тем более, когда это дает 50%-ное ускорение все программы.

Да, я еще упоминал, что вместо вектора, который засылается в функцию Multi() в качестве параметра у нас должна быть матрица, т.е. умножение немного усложнится. Пока можно тот вектор, который был раньше, расположить на диагонали новой матрицы.

    Поподробнее, пожалуйста, про матрицу в качестве параметра. Вы имеете в виду, что у параметра появится еще и мнимая часть (но тем не менее останется комплексным вектором) или все-таки станет буквально матрицей?
    Можно ли объяснить использование параметра в алгоритме. Дело в том, что входной вектор процедуры Multi(), по-видимому, является не сомножителем, а представляет собой "сырье", из которого там же приготавливается матрица-сомножитель, которая затем умножается на матрицу синусов/косинусов.
    Может быть, вы имеете в виду, что матрица-сомножитель будут формироваться в другом месте, а процедура будет получать ее уже в готовом виде в виде параметра?
     В любом случае хотелось бы, чтобы вы пояснили действия в процедуре Multi(), как в нынешнем ее виде, так и после изменения, на элементах матричной алгебры. Как выглядит этот алгоритм в матричной записи?

   Еще имею вопрос про выходные файлы. Отчего они называются как четные цифры? Почему пропущены нечетные? И что с ними потом делать?
   Чем больше цифра, тем сильнее мельчают интенсивности в файлах. На 5-7 кубитах файл 8.txt содержит одни нули. И даже на 8-кубитах содержит интенсивности лишь на уровне погрешности вычислений.
    Если строить график по "времени" во всем диапазоне (от 0 до 0.003), то в одном масштабе даже 6.txt не удастся увидеть – слишком мал.

[Pipa]

    С одним вопросом кажется сама разобралась - процедура Multi() выполняет преобразование:
U^T*D*U,
где D - диагональная, а U - унитарная
Унитарной матрицей служит таблица cosin, а диагональная представлена в параметре только своей диагональю.
    Видимо под матрицей в параметре вы имели в виду использование в качестве D матрицы любого типа, а не только диагональной.

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

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

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

Поподробнее, пожалуйста, про матрицу в качестве параметра. Вы имеете в виду, что у параметра появится еще и мнимая часть (но тем не менее останется комплексным вектором) или все-таки станет буквально матрицей?

Речь идет о той самой матрице ρ(0), о которой я писал в алгоритме. Можно ее взять вещественной, но лучше всего для общности сделать комплексной. Т.е. в качестве параметров в Multi() будут засылаться две матрицы (вещественная и мнимая часть ρ(0)), т.е. типа Multi(double **re_ro0, double **im_ro0).
Сама матрица ρ(0) должна задаваться в отдельной функции (или выбираться из стандартного набора, напр. равновесное состояние, cat-состояние, максимально-смешанное состояние и т.п.).

Можно ли объяснить использование параметра в алгоритме. Дело в том, что входной вектор процедуры Multi(), по-видимому, является не сомножителем, а представляет собой "сырье", из которого там же приготавливается матрица-сомножитель, которая затем умножается на матрицу синусов/косинусов.

Да нет, это не «сырье», а именно то, что умножается на синус/косинус:

Код:

for(i=0;i<N;i++)
         for(k=0;k<N;k++)
        {
         complex t(0.0,0.0);
         z1=c[i][k];
         z2=s[i][k];
         complex z(z1,-z2);
         t=z*plr[k];
         m2[i][k]=Re(t);
         rm1[i][k]=Im(t);   
        }

Здесь это вектор plr[k], который умножается на z, у нас вместо вектора должна быть матрица, причем лучше всего комплексная.

Для начала в defdiag.h можно создать отдельную функцию, которая будет формировать матрицу ro0. Например, Init_ro0(double *plr), создать в конструкторе новые матрицы для вещественной и мнимой части: double **re_ro0 и double **im_ro0, заполнить их нулями, а потом в новой функции их сформировать. В самом простом случае, чтобы отладить и получить совпадение с предыдущим вариантом, нужно оставить мнимую часть этой матрицы нулевой, а в вещественную матрицу заслать старый вектор на диагональные элементы матрицы re_ro0.

Код:

inline void BLOCK:: Init_ro0(double *plr)
{
int i;
for(i=0;i<N;i++)
re_ro0[i][i]= plr[i];
}

Затем в головной программе Main.cc перед основным циклом по времени вызвать эту функцию:
    bk->Init_ro0(bk->tm);

А тот кусок кода с умножением из функции Multi(), который я привел выше, нужно переписать, с учетом того, что вместо старого вектора plr[] будет уже комплексная матрица re_ro0[][], im_ro0[][]. Т.е. в матричной записи нужно прийти к тому варианту алгоритма, который я описывал в начале темы.

Еще имею вопрос про выходные файлы. Отчего они называются как четные цифры? Почему пропущены нечетные цифры? И что с ними потом делать?

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

Чем больше цифра, тем сильнее мельчают интенсивности в файлах. На 5-7 кубитах файл 8.txt содержит одни нули. И даже на 8-кубитах содержит интенсивности лишь на уровне погрешности вычислений.
    Если строить график по "времени" во всем диапазоне (от 0 до 0.003), то в одном масштабе даже 6.txt не удастся увидеть – слишком мал.

Да, номер максимальной когерентности не превосходит числа кубитов, например, в 3-кубитной системе будут только интенсивности 0 и 2-го порядков, в остальных файлах будут нули. И чем выше порядок когерентности, тем меньше она по абсолютной величине, это Вы правильно подметили. Но, повторюсь, нам эти интенсивности не нужны, мы будем считать квантовую запутанность и графики строить для нее.

[Pipa]

А тот кусок кода с умножением из функции Multi(), который я привел выше, нужно переписать, с учетом того, что вместо старого вектора plr[] будет уже комплексная матрица re_ro0[][], im_ro0[][]. Т.е. в матричной записи нужно прийти к тому варианту алгоритма, который я описывал в начале темы.

   Смотрите, что получилось.
После оптимизации кода старая процедура Multi() превратилась вот во что:

Код:

void BLOCK::MultiD(double *plr)
{
  double *re = new double[N];
  double *im = new double[N];

  for(int i=0; i<N; i++)
  { for(int k=0; k<N; k++)
    { re[k] = plr[k]*c[i][k];
      im[k] = -plr[k]*s[i][k];
    }

    for(int j=0; j<N; j++)
    { rm[i][j] = Scalar(re,c[j],N) - Scalar(im,s[j],N)
               + Scalar(re,s[j],N) + Scalar(im,c[j],N);
    }
  }
  delete[] im;
  delete[] re;
} 

Пояснения:
1) Функция Scalar(vec1,vec2,N) это не одноименный член класса Hami, а вызов ассемблерной функции вычисления скалярного произведения. Случайно совпали имена, но это не мешает работе, т.к. класс здесь другой.
2) Промежуточный продукт U^TD целиком не строится, а по мере надобности в том же цикле вычисляется рабочий столбец этого продукта. Благодаря чему, для работы дополнительно необходима только пара рабочих векторов для действительной и мнимой части рабочего столбца, а не полная матрица, как раньше.
3) Процедура выдает сумму действительной и мнимой части, но очевидно, что выкладывать результат можно по отдельности (вторую строку сумм не складывать, а укладывать в другой массив).
4) Функция Multi() переименована в MultiD(), подчеркивая в названии, что она работает с диагональю. Функцию, работающую с плотной комплексной матрицей, я назвала MultiC().

Аналогично, функция, работающая с комплексной матрицей, выглядит так:

Код:

void BLOCK::MultiC(double **re_ro0, double **im_ro0)
{
  double *re = new double[N];
  double *im = new double[N];
  
  for(int i=0; i<N; i++)
  { for(int k=0; k<N; k++)
    { re[k] = Scalar(re_ro0[i],c[i],N) - Scalar(im_ro0[i],s[i],N);
      im[k] = -Scalar(re_ro0[i],s[i],N) + Scalar(im_ro0[i],c[i],N);
    }

    for(int j=0; j<N; j++)
    { rm[i][j] = Scalar(re,c[j],N) - Scalar(im,s[j],N);
      rm1[i][j] = Scalar(re,s[j],N) + Scalar(im,c[j],N);
    }
  }
  delete[] im;
  delete[] re;
} 

Здесь действительная и мнимая части результата помещаются в матрицы rm[][] и rm1[][], сответственно.

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

Понял, спасибо, если результаты совпадают с исходными, то все отлично. С комплексной матрицей насколько дольше считается?
Кстати, Вы подсказали хорошую идею – сделать несколько вариантов функции Multi(), и в зависимости от типа начальной матрицы плотности выбирать тот или иной вариант, чтобы не считать каждый раз с комплексной матрицей. Нужно тогда еще третий вариант сделать для вещественной матрицы плотности ρ(0), напр. MultiR(double **re_ro0), он, похоже, чаще всего будет использоваться.

[Pipa]

Понял, спасибо, если результаты совпадают с исходными, то все отлично.

    Результаты с исходными совпадают полностью, за исключением 8.txt для 8-кубитной системы (при меньшем числе кубитов там одни нули). Значения интенсивностей там порядка е-34, из-за чего становятся видны ошибки округления (из-за ограниченности длины мантиссы теряются малые слагаемые на фоне большой суммы). У функции Scalar() погрешность меньше, т.к. складывет она в сумматор типа long double (лишних 11 двоичных разрядов, чем у типа double).

С комплексной матрицей насколько дольше считается?

    Почти вдвое дольше. Это очевидно из числа циклов. MultiC() вызывает Scalar() 8N² раз, а MultiD() только 4N² раза (обычными умножениями на этом фоне можно принебречь).

[Pipa]

    Еще один резерв скорости.
Матрица плотности, получающаяся в результате процедуры MultiC(), ведь всегда эрмитова? А, значит, фигли ее всю считать – достаточно верхнего треугольника! Модификация кода ничтожная – крутим внутренний цикл не до <N, а до <=i, а результат заталкиваем сразу в обе половинки матрицы (меняя местами только действительные и мнимые значения).
    Вроде бы, не теряя общности, получаем ускорение процедуры почти в два раза.
    Аналогичная Multi() задача оказалась и в процедуре cosin(), она тоже была оптимизирована тем же способом.

    Итого рекорд на сегодняшний день:

Кубиты      Время всей программы
     5                 0 sec
     6                 1 sec
     7                 8 sec
     8               47 sec (0.8 min)
     9             396 sec (6.6 min)

(Pentium-4 3.2 GHz, FSB 400 MHz).
расчет по MultiD().

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

Отлично! Так, глядишь, мы и 9 кубитов сможем осилить .

Значения интенсивностей там порядка е-34,

Это чистые нули, больше 20 знаков после запятой вообще нет смысла учитывать, число «пи» там только 20 знаков имеет.

Матрица плотности, получающаяся в результате процедуры MultiC(), ведь всегда эрмитова? А, значит, фигли ее всю считать – достаточно верхнего треугольника! Модификация кода ничтожная – крутим внутренний цикл не до <N, а до <=i, а результат заталкиваем сразу в обе половинки матрицы (меняя местами только действительные и мнимые значения).

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