Meneses Claudio
email: bmeneses@uch.edu.pe

Brian Alvarado-Diaz
email: witman.alvarado-diaz@ieee.orglimaperú

Witman Zamora
email: ezamora@uch.edu.pelima

Avid Eduardo Roman-Gonzalez
email: perúavid.roman-gonzalez@ieee.org

Implementation of A Wireless System for the Processing and Comparison of Cerebral Waves of Patients with Amyotrophic Lateral Sclerosis through Matlab Identifying Their Basic Needs

Keywords: Brainwaves, Amyotrophic Lateral Sclerosis (ALS), Wireless System

At present diverse illnesses exist between soft, mortal and silent, in which one some studies of the biomedical try to minimize and/or help to fight these illnesses. Due to this, in this paper presents a processing and comparison system of brainwaves, the patients with Amyotrophic Lateral Sclerosis, they are not able to move and/or use any of their body extremities so could communicate their basics needs by the brainwaves produced by their brains, taking apart the communication by the flickers, showing a more reliable and exact alternative. The following work consists of the design of a Wireless System able to show the brainwaves produced by the patient with Amyotrophic Lateral Sclerosis and also process and compare to get the meaning of each of them identifying the basics needs; moreover, it will benefit the family because it's going to exist a nonverbal communication with the patient.

I. INTRODUCTION

Today, medicine has evolved exponentially because of the great demand for studies and the needs of people, but also as presented in a positive development have also it has discovered new diseases between soft, mortal and loudness. The biomedical tries to contribute to the prevention, recovery, healing and in some cases in the implementation of sophisticated medical equipment to help society.

One of the mortal and silent illnesses is Amyotrophic Lateral Sclerosis, it's a degenerative illness that attacks mainly the nervous system, affecting the limbs of patients that suffer this illness, until the patient loses complete mobility of body and head, leaving only the movement of the eyes and eyelids.

The communication between the patients that suffer this illness with their relatives is very limited, the only existing method is by means of the chart of letters in which the only guide is the number patient's blinks in which it achieved known basic necessities such as: food, fatigue, acclimatization, etc., before this happens a great question of how to improve the communication of the patients with this illness.

Nowadays there are communication system for the people with ALS but most of them are using the eyes or the eyelids, it means that the patient has to use the movement of its eye erring because according some investigations we can keep opening our eyes just for 3 minutes, it depends where we are. Those communication system have many mistakes affecting mainly to the patient with ELA.

This thesis presents an additional method with the mayor efficient because it will be joined directly with the patients and also no works with the numbers of flickers, on the other hand, works with the brainwaves, with the aim of processing and comparing the brainwaves of patients with people who do not suffer this disease, to get the basic needs or needs at that time from the patient and be a safer and more efficient method for non-verbal communication with their relatives.

II. AMYOTROPHIC LATERAL SCLEROSIS

Amyotrophic Lateral Sclerosis is a degenerative disease of motor neurons in the spinal cord, brainstem or motor surface. This disease mainly affects the cells located in the lateral part of the body; this damage to cells are showed in muscle weakness. For this reason, patients' quality of life and their survival time can vary radically, so a medical procedure must be performed. Numerous researchers have focused on the field of treatment for neurodegenerative diseases, with the main objective of medical care towards these patients to improve their quality of life and thus obtain that survival time as long as possible. This disease is very influential in many countries due to its great impact of specializations based on neurodegenerative diseases. In order to improve the quality of life of these patients.

The next graphic shows the number or cases that happened in Peru. As there is in the graphic, Lima has the more incident of cases than other regions because most of the other regions the patients don't want to register theirselves. Also according to some rearchers there is a graphic that shows the years range where has more incident this desease. In other words between 30 to 59 years old, the incident of number of cases are higher than other years.

III. MATERIAL AND METHODS

To understand this article, I've structured the steps to get finished the general objective, including searching information,

A. Current Situation of the Communication Systems

Nowadays, there are several communication systems for patients with rare diseases or with dysfunctional motor. Many of these systems are used to improve the quality of life of patients with ALS, some of them are sophisticated electronic devices, they are effective but also very expensive alternative; in others types only required a table in which letters will be placed in an orderly manner and the patient by the movement of their eyes will indicate the direction to which it points; Therefore it needs a person who can write all the words and / or numbers that the patient with ALS points with the eyes. 1) Communication Type board ETRAN The ETRAN communication system consists in having a big board of a strong material because there will stick and / or adhere groups of words each one with different colors so that the patient with a motor dysfunction and / or brain dysfunctions, one of them known as ALS, using only the eyes can communicate by selecting the pictures and letters forming sentences. This communication system is really used in Peru because of its simple construction, in addition shows an evolutionary development [START_REF]ELA Mexico[END_REF] because the patient with ALS with the practice will be able to communicate much faster.

As can be seen in Figure 1 in this way the ETRAN board is constructed, successively that is positioned the listener person in charge of saying and / or confirming the color of the group and then the letter or number chosen by the patient with ALS.

The operation of the system is very simple, first the listener is placed next to the ETRAN board in front of the patient with ALS as shown in Figure 2.

Then through the movement of the eyes, the listener can recognize the group of letters or numbers selected by the patient, because the listener is staring at the patient face to face, can easily locate where the eyes are pointing.

When choosing the group of letters or numbers differentiated by the color, the listener will confirm what was established by the patient by repeating the color so that the patient can choose the letter or number, just as shown in Figure 2 each letter and number of each Group has colors, from which the patient will choose the color and thus can identify the letter or number that wanted to indicate. When the communication system has errors on the part of the patient, in some cases the same ETRAN table has a small section that says "erase", in other cases the patient only closes his eyes for an average of 4 seconds to indicate error [START_REF] Arboleda | Diseño y construcción de un prototipo de interfaz cerebro-computador para facilitar la comunicación de personas con discapacidad motora[END_REF][START_REF] Tseng V | Transparent: brain computer interface and social architecture[END_REF][START_REF] Feres | El método de las necesidades básicas insatisfechas (NBI) y sus aplicaciones en América Latina[END_REF].

2) Ocular communication

The ocular communication system is an optional language for people with motor disabilities, in other words, people who use the eyelids, eyes and eyebrows, being widely used by patients with Amyotrophic Lateral Sclerosis or other diseases with motor disabilities.

The ocular communication system was created by Marcelo Quintero, who at 37 years was diagnosed with ALS, thus thinking of an idea to improve his quality of life because this disease is degenerative. So, he thought of a communication system based on an ocular language using the eyelids, eyes and eyebrows, thus providing a non-invasive alternative of communication. The alternative ocular language was born from the idea of how an ocular movement could be translated to a sign that indicates a letter or number, thus allowing a communication similar to the verbal one. For this, 3 levels were prepared in which patients according to practice could improve their technician, therefore, communicate much faster. Below are 3 levels:

• First Level: Each eye movement will mean a number of which has a respective meaning. For example, following figure 3 belonging to the communication by the number of blinks; an ocular movement will mean a number, for this example eyes towards the right would have the number 2 position as meaning as shown in figure 4, having the meaning "I'm cold or I want my coat".

• Second level: At this level the eye communication system becomes more complex because it consists of combining 2 or more eye movements to form numbers greater than 2 to more numbers, fundamentally must have a table of actions and / or letters and / or numbers so that the patient with ALS and the person in his / her care can communicate. In case they do not have the table, with practice the patient with ALS and the person in charge of their care will be able to communicate, because the patient memorizes the contents of each number and the caregiver, each eye movement to what number does it belong

• Third level or last level: At this level it is necessary to learn the letters corresponding to the numbers formed by each eye movement by the patient with ALS, in order to be able to form words and sentences, in other words to be able to dialogue in a non-verbal with another person.

3) Communication by the number of blinks

The communication system by board of squares or also known by the number of flickers, as shown in figure 4, was used many times ago by many patients with Amyotrophic Lateral Sclerosis in Peru, this system had as purpose to use the number of flickers to be able to indicate the column and then the row, and that process obtained the letter or number and built the word letter by letter, taking a considerable time, taking into account that there were many failures due to the voluntary flicker that each person has by default. The patient with ALS could choose the letter or number that want to communicate and then confirmed by the same mode; in other cases to confirm only kept their eyes open and when the answer was contrary they closed the eyes for a period of time; so the listener or the assistant could know if it was right or wrong. [START_REF] Lee | A comparative study of wireless protocols: Bluetooth, UWB, ZigBee, and Wi-Fi[END_REF] Due to the time taken for each word, it was thought in a board for the same purpose but in this case they are actions, through the ease of images the patient with ALS was able to provide information about their needs basic or any action that desired in a required time. The wireless system will consist of various parameters for its operation; starting mainly with the frequency of brain waves in which we will identify the basic needs, then the device that will receive the brain waves in the same way will transmit in real time that data, taking into account its transmission power, finally ending with the receiving equipment being a portable computer, which will receive and process the data received.

1) Frequency of use of brain waves to identify basic needs:

The brainwaves are divided into 4 types according to their frequency level: Knowing that each type of brain wave has a specific meaning, it will be much easier to identify in which frequency range the patient with ALS is.

These measurements serve to be clear if the wireless system to be implemented will not have reception problems by NeuroSky Brainwave.

The brainwaves are formed by each movement or feeling that each person presents, being thus its representation of manner of signals; as shown in figure 6 are some examples of how the brain waves are drawn according to each action of the person

2) Standard Dimensions of a Patient Unit

The patient unit is the conglomerate formed by the space of the room, the furniture conformed by the bed, a mattress, pillow with cover, lateral handrails so that the patient does not fall and a crank to be able to regulate the height of the bed, and finally all the furniture equipment that is observed in the figure 7 and the material and equipment are provided to the patient for his daily activities like the cleaning, feeding, etc .; the patient uses it for his / her stay in the hospital, it means that the number of patient units must be equal to the number of beds within that room. As shown in Figure 7, a patient unit is shown next to the respective elements that it must has.

It is common for hospitals to separate rooms depending on the type of illness, age, type of social class and number of beds, therefore the organization and arrangement of the patient unit will depend on the hospital.

The dimensions of separation of the beds of each patient is from 1 to 1.20 meters, means that for the use of the wireless system will not have problems because it will be close to the patient, without obstructions by means.

C. Steps to Follow for the Implementation of the Wireless System

1) Connecting the NeuroSky Brainwave to the laptop The NeuroSky Brainwave, when turned on, can be denoted by a blue LED as shown in Figure 8 and emits its own Bluetooth network so that other devices having equally bluetooth can be paired with the NeuroSky Brainwave. When we observe that the blue LED lights correctly on the NueroSky Brainwave we proceed to pair from the laptop using its Bluetooth. When we add in a Bluetooth device we will see 2 option with the name of MindWave Mobile both of the same name but based on the manual that comes with the device just have to pair with the one with the icon the headphones.

When we have already paired the sign of the headphones, it means we can transmit brain waves to the laptop.

Note: On some laptop at the time of pairing requests a key, according to the manual, the default password is 0000.

2) Receiving Data from NeuroSky Brainwave Software

The NeuroSky Brainwave does not have any element that shows that it is transmitting data; the only sign that shows that it is on and supposedly transmitting is its blue led. Because of this we will test it with the laptop, the NeuroSky Brainwave comes with a CD of its software, and it has applications to know if it is transmitting data.

When we install the program, it will show a window, in which the icon of the installation and a window with 2 applications: BrainWave Visualizer and SpeedMath.

After installing the software we will proceed in the correct positioning of the NeuroSky Brainwave, as shown in figure 9 is as it will be seen from front, as it is observed the frontal sensor detects the electrical signals of the brain known as brain waves, it is possible to emphasize that the sensor Frontal must be interacting with the skin, if in case it is not that will present errors in wanting to capture the brainwaves.

It should be clarified that the arm where the front sensor is located must be lowered 90 ° and thus be positioned.

In the same way as shown in figure 10, in this case it is in profile, where we realize that the EarClip is also connected in the inner part of the ear, this device will act as a ground and will help to filter all the Electrical noise of the body and environment so as to have a purer brain wave without distortions. Also shown is the opening where an AAA battery will be inserted to power the NeuroSky Brainwave.

When we have correctly located the NeuroSky Brainwave, we will proceed to open the application called BrainWave Visualizer to be able to visualize our brain waves on the computer screen as shown in figure 10, where in the upper right you can see the meter of signal; when it really means the NeuroSky Brainwave is connected and transmitting correctly. In the lower right we see 2 levels that indicate in real time how much attention and / or how much meditation we are experiencing with our brain. Then a sphere is displayed where the brain waves will be shown according to the brain interaction that is receiving the NeuroSky Brainwave. Finally the bars indicate the frequency of our brainwaves.

3) Connecting the NeuroSky Brainwave to MatLab

Knowing that MatLab is a very used and powerful software in terms of signal processing, we will connect the NeuroSky Brainwave to MatLab by means of commands that indicates the reception of the brainwaves and these can be drawn and visualized by the same software.

To connect in NeuroSky Brainwave to MatLab after having paired to the laptop we must know in which COM port is connected. It shows the port on which the laptop is receiving data. Note two COM ports that both carry the name of MindWave Mobile, we must choose the COM port that has an INCOMING address because we will receive the information. In my case the NeuroSky Brainwave is in the COM3 port, this information will be added in the configuration in the MatLab because the data reception input must be declared.

By executing the MatLab programming, previously using the NeuroSky Brainwave properly, we will observe the data in numbers mode that are received by the laptop and then a graph that is drawn according to the data that is being added. It can be seen in detail in Figure 12, how the brain waves are constructed.

IV. RESULTS

In the study of the current situation, I found that the current communication systems of patients with ALS are very inaccurate and deplorable in their data because it requires the physical effort of the patient such as flickers, eye movement, thus erring the information that the patient relief to transmit to his relatives.

In the study of the parameters for the operation of the system, I was able to identify which are the fundamental parts in the system, because each parameter of the equipment mainly affects the system. Knowing the dimensions of the patient's room, brainwave and laptop metrics, NeuroSky Brainwave information, optimal system performance can be achieved. Unfortunately, I can´t analyze a patient with ELA because they are very restricted in their visit, they just have visits by their relatives.

In the study steps for the implementation of the wireless system, I found the connection of the NeuroSky Brainwave to the laptop and then use MatLab software to draw the brainwave. Also the correct operation and position of the device so that the wireless system has no problem. As a result I got the steps to declare the NeuroSky Brainwave in MatLab.

V. DISCUSSION

The present paper confirms the communication to identify the basic needs of patients with ALS through the brain waves are one of the most efficient and modern communication alternatives for patients with ALS due to their degeneration of motor neurons that finally removing the movement of any part of the body, only with the ability to move the eyes and eyelids.

There are currently many jobs related to improving the quality of life of people suffering degenerative diseases due to their social impact attributed to them, the majority of these communication systems are using the eyes and movement of the same, providing no communication -verbal; also having possible errors that occur. The communication systems using the eyes imply that there must be an adequate environment also a prior knowledge about what communication will consist because there are several types, the most used is by means of boards of which the patient using the flickers will indicate the letter or numbers or actions. These communication systems improve the quality of life of patients with ALS but are not very effectives, causing dissatisfaction and stress to the patient with ALS who wants to indicate something. Many of the studies students to patients with ALS only to know where the disease comes from but not many are looking for alternative means to improve their quality of life.

In conclusion, in many cases, communication through the brain waves is an innovative technological system where studies are already being made based on these systems because of their great effectiveness in improving the quality of life of patients with degenerative diseases.

CONCLUSIONS

To concluded that currently the communication systems for patients in ALS are those that use eye movement and general use of eyes, eyes and eyebrows. Being very easy to use and create; what is required in conclusion is to improve the quality of life of patients with ALS.

To concluded that having the necessary parameters taking into account the brain waves, the patient unit, the device that will receive and send the cerebrals waves and the reception sensitivity of the receiving equipment, we note that they are key terms for the operation of the system, so that at the time of the final test there are no connection problems or terminologies.

To conclude that knowing the correct connection of the NeuroSky Brainwave, and also how to test it in a programming in the MatLab, we can get the brainwaves without any inconvenience. Also keep in mind that the NeuroSky Brainwave to work in MatLab software, should download the libraries that the same computer offers on its own website.

Figure 1 :Figure 1 :

 11 Figure 1:Table of Communication ETRAN (Numbers Included)

Figure 2 :Figure 4 :Figure 5 :

 245 Figure 2: Each ocular movement means a number

Figure 6 :

 6 Figure 6: Representations by brainwaves of for each action of the person

Figure 7 :

 7 Figure 7: Room of a Patient Unit

Figure 8 :

 8 Figure 8: Blue LED light showing that the device is working

Figure 9 :Figure 11 :Figure 10 :

 91110 Figure 9: Using the NeuroSky Brainwave Front View

Figure 12 :Figure 13 :

 1213 Figure 12: Construction of the brainwave in MatLab

TABLE 1 .

 1 DESCRIPTION OF THE BRAINWAVES

	Brainwaves Frequency	Electrical Level Function
	Delta	0 -4 Hz	10 to 50 μV	deep sleep
	Theta	4 -8 Hz	50 to 100 μV	meditation
	Alpha	8-12 Hz	100 to 150 μV	daydreaming


Qingxiang Feng

Senior Member, IEEE Yicong Zhou
email: yicongzhou@umac.mo

Color Face Recognition using High-Dimension Quaternion-based Adaptive Representation

Keywords: Sparse Representation, Collaborative Representation, Quaternion, Color Face Recognition

Recently, quaternion collaborative representationbased classification (QCRC) and quaternion sparse representation-based classification (QSRC) have been proposed for color face recognition. They can obtain correlation information among different color channels. However, their performance is unstable in different conditions. For example, QSRC performs better than than QCRC on some situations but worse on other situations. To benefit from quaternion-based e2-norm minimization in QCRC and quaternion-based e1-norm minimization in QSRC, we propose the quaternion-based adaptive representation (QAR) that uses a quaternion-based ep-norm minimization (1 ≤ p ≤ 2) for color face recognition. To obtain the high dimension correlation information among different color channels, we further propose the high-dimension quaternion-based adaptive representation (HD-QAR). The experimental results demonstrate that the proposed QAR and HD-QAR achieve better recognition rates than QCRC, QSRC and several state-of-the-art methods.

I. INTRODUCTION

F ACE recognition systems are relied on classification methods. Recently, representation-based classification methods [START_REF] Wright | Sparse representation for computer vision and pattern recognition[END_REF]- [START_REF] Feng | Two classifiers based on nearest feature plane for recognition[END_REF] have attracted more attention and obtained good performance in face recognition. Two typical representation-based classification methods are sparse representation-based classification (SRC) [START_REF] Wright | Sparse representation for computer vision and pattern recognition[END_REF] [START_REF] Wright | Robust face recognition via sparse representation[END_REF] and collaborative representation based classification (CRC) [START_REF] Zhang | Sparse representation or collaborative representation: Which helps face recognition?[END_REF]. For color face recognition, SRC and CRC independently use color channels of color face images. As a result, they fails to consider the structural correlation information among different color channels. The previous research works [START_REF] Pei | Efficient implementation of quaternion fourier transform, convolution, and correlation by 2-d complex fft[END_REF]- [START_REF] Barthélemy | Color sparse representations for image processing: review, models, and prospects[END_REF] proved that quaternion can obtain the information among different channels of the color images by representing three color channels as three imaginary parts. For example, [START_REF] Kolaman | Quaternion structural similarity: a new quality index for color images[END_REF] proposed a quaternion-based structural similarity index and used it to assess the quality of color images. [START_REF] Xu | Vector sparse representation of color image using quaternion matrix analysis[END_REF] proposed a vector sparse representation using quaternion matrix analysis and successfully applied it to various color image processing. [START_REF] Zou | Quaternion collaborative and sparse representation with application to color face recognition[END_REF] proposed two methods to obtain correlation information among red, green and blue channels for color face recognition. These two methods are named quaternion collaborative representation-based classification (QCRC) and

• To benefit both from quaternion-based e 2 -norm minimization in QCRC and quaternion-based e 1 -norm minimization in QSRC, we propose quaternion-based adaptive representation (QAR) method. QAR uses a quaternionbased e p -norm minimization (1 ≤ p ≤ 2) that can adaptively combine quaternion-based e 2 -norm minimization and quaternion-based e 1 -norm minimization. • To our best knowledge, the high-dimension quaternion has not been applied. Thus, we propose the highdimension quaternion-based adaptive representation (HD-QAR). It is able to obtain high-dimension correlation information among different color channels. • Because the product of high-dimension quaternion is noncommunicative, it is quite difficult to solve the highdimension quaternion-based optimization problem in HD-QAR. Thus, we propose the quaternion-based kernel matrix and two high-dimension-quaternion operators to solve optimization problem in HD-QAR be easy.

• Extensive experiments on six well-known color face databases show that the proposed methods have better performance than several state-of-the-art methods.

II. RELATED WORK

This section reviews the related works including quaternion, kernel and representation-based methods.

A. Quaternion Algebra

Hamilton proposed the quaternion space that has three imagery units i, j and k [START_REF] Hamilton | Ii. on quaternions; or on a new system of imaginaries in algebra[END_REF]. These three imagery units are defined as

i 2 = j 2 = k 2 = ijk = -1 (1)
A quaternion q is described by

q = q 0 + q 1 i + q 2 j + q 3 k (2)
where q 0 , q 1 , q 2 , q 3 are real numbers. The conjugate quaternion of q is

¯q = q 0 -q 1 i -q 2 j -q 3 k (3)
The modulus of quaternion q is

| q| = q ¯q = q 2 0 + q 2 1 + q 2 2 + q 2 3 ( 4
)
The addition operation of two quaternions ṗ and q is q + ṗ = (q 0 + p 0 ) + (q 1 + p 1 )i + (q 2 + p 2 )j + (q 3 + p 3 )k (5)

The multiplication operation of two quaternions ṗ and q is q ṗ = (q 0 p 0 -q 1 p 1 -q 2 p 2 -q 3 p 3 ) +(q 0 p 1 + q 1 p 0 + q 2 p 3 -q 3 p 2 )i +(q 0 p 2 -q 1 p 3 + q 2 p 0 + q 3 p 1 )j +(q 0 p 3 + q 1 p 2 -q 2 p 3 + q 3 p 0 )k

B. Kernel trick

The kernel techniques can map the data to a high-dimension space. Using the kernel trick, we can obtain the high dimension information without computing the mapping explicitly. In this paper, we use a well-known kernel called the Gaussian radial basis function (RBF) kernel. It can be represented as

k(a, b) = φ(a) T φ(b) = exp( -||a -b|| 2 δ ) (7)
where a and b are any two original samples, δ is a parameter.

In the kernel trick, φ( * ) is unknown. We can access the feature space only via k( * , * ) .

They are two typical sparse-based methods. A lot modifications of SRC/CRC have been proposed for various visual recognition tasks [START_REF] Xu | A two-phase test sample sparse representation method for use with face recognition[END_REF]- [START_REF]Kernel combined sparse representation for disease recognition[END_REF]. For example, Gao et al. denoted that the nonlinear high dimensional features are useful for sparse representation, and then proposed kernel-based sparse representation classification (KSRC) [START_REF] Gao | Kernel sparse representation for image classification and face recognition[END_REF], [START_REF] Gao | Sparse representation with kernels[END_REF],

min ||α|| 1 s.t. Φ(x) = Φ(X)α (10)
Wang et al. proposed kernel-based collaborative representation classification (KCRC) [START_REF] Wang | Kernel collaborative representation-based classifier for face recognition[END_REF], which is described as

min ||α|| 2 s.t. Φ(x) = Φ(X)α (11)
Zou et al. denoted that quaternion can help sparse representation to obtain the correlation information among different color channels. Then, they proposed quaternion sparse representation classification (QSRC) and quaternion collaborative representation based classification (QCRC). QSRC is described as

min || α|| 1 s.t. ẋ = Ẋ α (12)
QCRC is described as

min || α|| 2 s.t. ẋ = Ẋ α (13)

III. PROPOSED QAR

In this section, we propose a new method, called quaternionbased adaptive representation (QAR). We first introduce its motivation. Then, we define and analyze its objective function. Afterwards, the solution of this objective function is given. Last, we describe the classification rule and summary this method.

A. Motivation

To consider the structural correlation information among different color channels, [START_REF] Zou | Quaternion collaborative and sparse representation with application to color face recognition[END_REF] proposed two methods for color face recognition. They are quaternion collaborative representation-based classification (QCRC) using e 1 -norm minimization and quaternion sparse representation-based classification (QSRC) using e 2 -norm minimization. Table II lists some results copied from [START_REF] Zou | Quaternion collaborative and sparse representation with application to color face recognition[END_REF]. From Table II and [START_REF] Zou | Quaternion collaborative and sparse representation with application to color face recognition[END_REF], we find that QSRC obtains the better performance than QCRC on some situations while QCRC obtains the better performance than QSRC on other situations. This means that quaternion-based e 2 -norm minimization and quaternion-based e 1 -norm minimization both have their own advantages. Motivated by this, this paper designs a quaternion-based e p -norm minimization (1 ≤ p ≤ 2) to solve the sparse optimization problem such that it can benefit from quaternion-based e 2 -norm minimization in QCRC and quaternion-based e 1 -norm minimization in QSRC. We give an example in Fig. 1 to describe the geometric interpretation of e p -norm, e 2 -norm and e 1 -norm.

B. Objective function of QAR

V = V 0 + V 1 i + V 2 j + V 3 k . ( V ) is equal to ( V ) :=     V 0 -V 1 -V 2 -V 3 V 1 V 0 -V 3 V 2 V 2 V 3 V 0 -V 1 V 3 -V 2 V 1 V 0     (16) Given a vector-based quaternion v = v 0 + v 1 i + v 2 j + v 3 k. ℵ( v) is equal to ℵ ( v) = v T 0 v T 1 v T 2 v T 3 T ( 17
)
Using the two quaternion-based operations ( * ) and ℵ( * ), the objective function in Eq. ( 14) can be rewritten as

min ||ℵ( α)|| p s.t. ℵ( ẋ) = ( Ẋ)ℵ( α) (18)
and the objective function in Eq. ( 15) can be rewritten as

min || ( Ẋ)Diag(ℵ( α))|| * s.t. ℵ( ẋ) = ( Ẋ)ℵ( α) (19)
Next, we will explain why Eq. ( 19) can be treated as e p -norm (1 ≤ p ≤ 2).

C. Analysis of objective function

This section analyzes the objective function and explain why min || ( Ẋ)Diag(ℵ( α))|| * in Eq. ( 19) is equivalent to min ||ℵ( α)|| p in Eq [START_REF] Gao | Kernel sparse representation for image classification and face recognition[END_REF]. In order to explain the proposed objective function can be treated as e p -norm minimization (1 ≤ p ≤ 2), we give two examples. Suppose that the subjects (columns) of ( Ẋ) are orthogonal and different from each other, that is, ( Ẋ) T ( Ẋ) = I , where I is an identity matrix. We define Ω = || ( Ẋ)Diag(ℵ( α))|| * and its decomposition is 20), we know that the objective function in Eq. ( 19) is equivalent to e 1 -norm minimization problem. On the contrary, we suppose that the subjects (columns) of ( Ẋ) are same to the first column of ( Ẋ) , that is ( Ẋ) T ( Ẋ) = 11 T ( 1 is a vector that all elements are one). We define the first column of ( Ẋ) as ( Ẋ) 1 and The decomposition is

Ω = T r ( Ẋ)Diag(ℵ( α)) T ( Ẋ)Diag(ℵ( α)) 1 / 2 = T r (Diag(ℵ( α))) T (Diag(ℵ( α))) 1 / 2 = ||ℵ( α)|| 1 (20) Based on Eq. (
Ω = || ( Ẋ) 1 1 T Diag(ℵ( α))|| * (21)
Suppose that the rank of a matrix A ∈ R m×n is r. We can obtain the inequalities [START_REF] Golub | Matrix computations[END_REF] [28] as follows.1

||A|| F ≤ ||A|| * ≤ √ r||A|| F (22)
where || • || F is the Frobenius norm. Because the rank of ( Ẋ)

1 1 T Diag(ℵ( α)) is one (r = 1), Ω = || ( Ẋ)Diag(ℵ( α))|| * can be re-wrote as Ω = || ( Ẋ) 1 1 T Diag(ℵ( α))|| * = || ( Ẋ) 1 1 T Diag(ℵ( α))|| F = || ( Ẋ) 1 || 2 ||1 T || 2 ||Diag(ℵ( α))|| 2 = ||ℵ( α)|| 2 (23)
In the real applications, the subjects (columns) of data set ( Ẋ) are not different each other or same to each other. Consider the Eqs. ( 20) and ( 23), our objective function could be treated as a combination of the quaternion-based e 1norm minimization in QSRC and quaternion-based e 2 -norm minimization in QCRC. That is, Eq. ( 19) can be treated as quaternion-based e p -norm. Our objective function could benefit from both quaternion-based e 1 -norm minimization and quaternion-based e 2 -norm minimization based on the correction of the subjects (columns) in the data set ( Ẋ) .

D. Solution of objective function

Based on the optimization methods in [START_REF] Candès | Robust principal component analysis?[END_REF] and [START_REF] Liu | Robust subspace segmentation by low-rank representation[END_REF], we use the inexact augmented Lagrange multipliers (IALM) [40] method to solve the proposed e p -norm minimization. We convert the minimization problem in Eq.( 19) to the following optimization problem as min Z,z,ℵ( α)

||Z|| * + ||z|| 1 s.t. z = ℵ( ẋ) -( Ẋ)ℵ( α) Z = ( Ẋ)Diag(ℵ( α)) (24)
We further convert the optimization problem in Eq. ( 24) to a augmented Lagrange multiplier problem as follows.

L(Z, z, ℵ( α)) = λ||Z|| * + m T 1 (ℵ( ẋ) -( Ẋ)ℵ( α) -z) + ||z|| 1 + T r[M T 2 (Z -( Ẋ)Diag(ℵ( α)))] + u 2 (||a -( Ẋ)ℵ( α) -z|| 2 2 + ||Z -( Ẋ)Diag(ℵ( α))|| 2 F ) (25)
where u > 0 is a parameter, m 1 and M 2 are Lagrange multipliers. We can alternatively optimize the variates Z, z and ℵ( α) by fixing others. The detailed optimization processes of Z, z and ℵ( α) are described as follows. Update Z by fixing z and ℵ( α) . It is to solve the following minimization problem as

Z * = arg min Z L(Z, z, ℵ( α)) = arg min Z λ||Z|| * + T r(M T 2 Z) + u 2 ||Z -( Ẋ)Diag(ℵ( α))|| 2 F = arg min Z λ u ||Z|| * + 1 2 ||Z -( ( Ẋ)Diag(ℵ( α)) -1 u M 2 )|| 2 F ( 26
)
We can approximately solve the above minimization problem in Eq. ( 26) using the singular value thresholding (SVT) operator [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF]. Then, update ℵ( α) by fixing z and Z. It is to solve the following minimization problem as.

ℵ( α) * = arg min

ℵ( α) L(Z, z, ℵ( α)) = arg min ℵ( α) -m T 1 Jℵ( α) -T r(M T 2 JDiag(ℵ( α))) + u 2 (ℵ( α) T J T Jℵ( α) -2(ℵ( ẋ) -z) T Jℵ( α) -2T r(Z T JDiag(ℵ( α))) + T r((JDiag(ℵ( α))) T JDiag(ℵ( α)))) = arg min ℵ( α) u 2 ℵ( α) T (J T J + Diag(diag(J T J)))ℵ( α)- (J T m 1 + uJ T J(ℵ( ẋ) -z) + diag(M T 2 J + uZ T J)) T ℵ( α) (27)
where J = ( Ẋ). The above minimization problem in Eq. ( 27) can be solved by

ℵ( α) * = (J T J + Diag(diag(J T J))) -1 J T ( 1 u m 1 + J(ℵ( α) -z)) +(J T J + Diag(diag(J T J))) -1 diag(J T ( 1 u M 2 + Z)) (28)
Next, update z by fixing Z and ℵ( α). It is to solve the following minimization problem as

z * = arg min z L(Z, z, ℵ( α)) = arg min z ||z|| 1 -m T 1 z + u 2 ||a -( Ẋ)ℵ( α) -z|| 2 F = arg min z 1 u ||z|| 1 + 1 2 ||z -(a -( Ẋ)ℵ( α) + 1 u m 1 )|| 2 2 (29)
We can solve the minimization problem in Eq. ( 29) by using the soft thresholding (shrinkage) operator [START_REF] Hale | Fixed-point continuation for \ell 1-minimization: Methodology and convergence[END_REF]. After updating Z, z and ℵ( α), the multipliers m 1 , M 2 can be updated by

m 1 = m 1 + u(ℵ( ẋ) -( Ẋ)ℵ( α) -z) M 2 = M 2 + u(Z -( Ẋ)Diag(ℵ( α))) (30)
We update the parameter u by u = min(ρu, u max ), where ρ, u max are two constant values.

E. Classification

After obtaining the results of the optimization problem of the objective function, we calculate the distance between testing sample x and each class by

d c = ||ℵ( ẋ) -( Ẋc )ℵ( αc )|| (31)
Finally, the testing sample will be classified into the class with the minimization distance by

c * = arg min(d c ) (32)
The detailed classification procedures of QAR are summarized in Algorithm 1.

IV. PROPOSED HD-QAR

This section proposes high-dimension quaternion based adaptive representation (HD-QAR). We first introduce the motivation and its objective function. The solution of this objective function and the classification rule are then presented. Finally, we compare the proposed methods with several stateof-the-art methods.

L(Z, z, ℵ( α)) = λ||Z|| * + ||z|| 1 + m T 1 (a -( Ẋ)ℵ( α) -z) + T r[M T 2 (Z -( Ẋ)Diag(ℵ( α)))] + u 2 (||ℵ( ẋ) -( Ẋ)ℵ( α) -z|| 2 2 + ||Z -( Ẋ)Diag(ℵ( α))|| 2 F ) 3: while until convergence do 4:
Update Z by fixing the others as Z * = arg min

Z λ u ||Z|| * + 1 2 ||Z -( ( Ẋ)Diag(ℵ( α)) -1 u M 2 )|| 2 F 5:
Update ℵ( α) by fixing the others as

ℵ( α) * = ( ( Ẋ) T ( Ẋ) + Diag(diag( ( Ẋ) T ( Ẋ)))) -1 ( Ẋ) T ( 1 u m 1 + ( Ẋ)(ℵ( ẋ) -z)) +( ( Ẋ) T ( Ẋ) + Diag(diag( ( Ẋ) T ( Ẋ)))) -1 diag( ( Ẋ) T ( 1 u M 2 + Z)) 6:
Update z by fixing the others as z * = arg min

z 1 u ||z|| 1 + 1 2 ||z -(ℵ( ẋ) -( Ẋ)ℵ( α) + 1 u m 1 )|| 2 7:
The

A. Motivation and Objective Function

To better describe the motivation, we give an example of color face images shown in Fig. 2. A color face image has red, green and blue channels. The quaternion can obtain structural correlation information among different color channels. [START_REF] Gao | Kernel sparse representation for image classification and face recognition[END_REF], [START_REF] Gao | Sparse representation with kernels[END_REF], [START_REF] Zhang | Kernel sparse representation-based classifier[END_REF] presented high dimension information in kernel space are helpful for sparse. Thus, we want to design a kernel quaternion transformation that can obtain the high dimension correlation information among different color channels. Suppose that there is a mapping function, Φ(.) : R q → R Q (q << Q) , which maps a vector-based quaternion ẋ and a matrix-based quaternion Ẋ to the high dimensional feature space as x → Φ( ẋ) X → Φ( Ẋ) . Using this mapping function, the objective function of HD-QAR is as follows.

min || α|| p s.t. Φ( ẋ) = Φ( Ẋ) α ( 33
)
where 1 ≤ p ≤ 2.

B. Two kernel-based quaternion operations

Given a matrix-based quaternion

V = V 0 + V 1 i + V 2 j + V 3 k , the high dimension matrix quaternion is φ( V ) = φ(V 0 ) + φ(V 1 )i + φ(V 2 )j + φ(V 3 )k . Given a vector-based quaternion v = v 0 +v 1 i+v 2 j+v 3 k , the high dimension matrix quaternion is φ( v) = φ(v 0 ) + φ(v 1 )i + φ(v 2 )j + φ(v 3 )k . Because φ( V )
and φ( v) contain the complex number, we need to design two operations to transform φ( V ) and φ( v) to the real number such that the optimization problem in Eq.( 33) can be computed easily. These operations are

φ ( V ) :=     φ(V 0 ) -φ(V 1 ) -φ(V 2 ) -φ(V 3 ) φ(V 1 ) φ(V 0 ) -φ(V 3 ) φ(V 2 ) φ(V 2 ) φ(V 3 ) φ(V 0 ) -φ(V 1 ) φ(V 3 ) -φ(V 2 ) φ(V 1 ) φ(V 0 )     (34)
and

ℵ φ ( v) = φ(v T 0 ) φ(v T 1 ) φ(v T 2 ) φ(v T 3 ) T (35)

C. Solution of objective function

Using the proposed two kernel-based quaternion operations, the proposed objective function of HD-QAR can be represent-Algorithm 2 Proposed High-Dimension Quaternion-based Adaptive Representation (HD-QAR) Method Require: The original training set X with L samples and a testing sample x. Initialize the value of S, s, α, m 1 , M 2 , u, ρ, ε and u max Ensure: Label of the testing sample x.

1: Use X and x to constitute quaternion-based kernel matrix and kernel vector as

K =   k ℵ φ ( ẋ1 ), ℵ φ ( ẋ1 ) ... k ℵ φ ( ẋ1 ), ℵ φ ( ẋL ) ... ... ... k ℵ φ ( ẋL ), ℵ φ ( ẋ1 ) ... k ℵ φ ( ẋL ), ℵ φ ( ẋL )   k = k ℵ φ ( ẋ1 ), ℵ φ ( ẋ) ... k ℵ φ ( ẋL ), ℵ φ ( ẋ)
L(S, s, γ) = λ||S|| * + ||s|| 1 + y T 1 (k -Kγ -s) + T r[Y T 2 (S -KDiag(γ))] + u 2 (||k -Kγ -s|| 2 2 + ||S -KDiag(γ)|| 2 F )
4: while until convergence do Update S by fixing the others with S * = arg min

S λ u ||S|| * + 1 2 ||S -(KDiag(γ) -1 u M 2 )|| 2 F 6:
Update γ by fixing the others with

γ * = (K T K + Diag(diag(K T K))) -1 K T ( 1 u m 1 + K(k -s)) + (K K + Diag(diag(K T K))) -1 diag(K T ( 1 u M 2 + S)) 7:
Update s by fixing the others with s * = arg min

s 1 u ||s|| 1 + 1 2 ||s -(x -Kγ + 1 u m 1 )|| 2 8:
The multipliers m 1 , M 2 are updated by

m 1 = m 1 + u(k -Kγ -s) and M 2 = M 2 + u(S -KDiag(γ)) 9:
Update the parameter u by u = min(ρu, u max )

Because the product of high-dimension quaternion are noncommunicative, solving the high-dimension quaternion-based optimization problem in HD-QAR is quite challenging. Thus, we propose quaternion-based kernel matrix K and quaternionbased kernel vector k as follows. Suppose that K = φ ( Ẋ) T φ ( Ẋ) and k = φ ( Ẋ) T ℵ φ ( ẋ), which can be computed as follows:

K = φ ( Ẋ) T φ ( Ẋ) =   k ℵ φ ( ẋ1 ), ℵ φ ( ẋ1 ) ... k ℵ φ ( ẋ1 ), ℵ φ ( ẋL ) ... ... ... k ℵ φ ( ẋL ), ℵ φ ( ẋ1 ) ... k ℵ φ ( ẋL ), ℵ φ ( ẋL )   (37) and k = φ ( Ẋ) T ℵ φ ( ẋ) = k ℵ φ ( ẋ1 ), ℵ φ ( ẋ) ... k ℵ φ ( ẋL ), ℵ φ ( ẋ) (38)
Using the Eqs.( 37) and ( 38), the Eq.( 36) is rewritten as min S,s,γ

||S|| * + ||s|| 1 s.t. s = k-Kγ S = KDiag(γ) (39)
The minimization problem in Eq.( 39) is then transformed into the following augmented Lagrange multiplier problem as

L(S, s, γ) = λ||S|| * + ||s|| 1 + y T 1 (k -Kγ -s) + T r[Y T 2 (S -KDiag(γ))] + u 2 (||k -Kγ -s|| 2 2 + ||S -KDiag(γ)|| 2 F ) (40
) where u > 0 is a parameter, m 1 and M 2 are Lagrange multipliers. Similar to the proposed QAR, variables S, s and γ in Eq.( 40) can be optimized alternatively when the other two variables are fixed. Update S when s and γ are fixed, which is equivalent to solve the minimization problem

S * = arg min E L(S, s, γ) = arg min S λ||S|| * + T r(M T 2 E) + u 2 ||S -KDiag(γ)|| 2 F = arg min S λ u ||S|| * + 1 2 ||S -(KDiag(γ) -1 u M 2 )|| 2 F (41)
We use the singular value thresholding (SVT) operator [START_REF] Cai | A singular value thresholding algorithm for matrix completion[END_REF] to approximately solve the minimization problem in Eq.( 41). Afterwards, update γ when s and S are fixed, which can be

= arg min γ -m T 1 Kγ -T r(M T 2 KDiag(γ)) + u 2 (γ T K T Kγ -2(k -s) T Kγ -2T r(S T KDiag(γ)) +T r((KDiag(γ)) T KDiag(γ))) = arg min γ u 2 γ T (K T K + Diag(diag(K T K)))γ- (K T m 1 + uK T K(k -s) + diag(M T 2 K + uS T K)) T γ (42)
The above minimization problem can be easily solved by

γ * = (K T K + Diag(diag(K T K))) -1 K T ( 1 u m 1 + K(k -s)) +(K T K + Diag(diag(K T K))) -1 diag(K T ( 1 u M 2 + S)) (43
) Next, update s when S and γ are fixed, which can be solved as follows.

s * = arg min s L(S, s, γ) = arg min s ||s|| 1 -m T 1 s + u 2 ||k -Kγ -s|| 2 F = arg min s 1 u ||s|| 1 + 1 2 ||s -(k -Kγ + 1 u m 1 )|| 2 2 (44)
We use the soft thresholding (shrinkage) operator [START_REF] Hale | Fixed-point continuation for \ell 1-minimization: Methodology and convergence[END_REF] to solve the solution of the minimization problem in Eq.( 44). After updating S, s and γ, the multipliers m 1 , M 2 are updated by

m 1 = m 1 + u(k -Kγ -s) M 2 = M 2 + u(S -KDiag(γ)) (45)
Then, update u by u = min(ρu, u max ), where ρ, u max are two constants.

D. Classification

Using the sparse coefficient γ, we compute the distance between testing sample and the c th class as

d c (x) = ||k -K c α c || ||α c || (46)
The rule of HD-QAR is to select the class with the minimum one as

min c * d c (x) , c = 1, 2, ..., M (47)
The complete classification procedure of HD-QAR is shown in Algorithm 2.

E. Proposed methods vs Well-known methods

The main difference between the proposed methods and the well-known methods is the way how to obtain the sparse coefficients. Table III summarizes and compares the proposed QAR and HD-QAR with these methods.. SRC [START_REF] Wright | Sparse representation for computer vision and pattern recognition[END_REF] solves e 1 -norm minimization problem. CRC [START_REF] Zhang | Sparse representation or collaborative representation: Which helps face recognition?[END_REF] solves e 2norm minimization problem. KSRC [START_REF] Gao | Kernel sparse representation for image classification and face recognition[END_REF], [START_REF] Gao | Sparse representation with kernels[END_REF], [START_REF] Zhang | Kernel sparse representation-based classifier[END_REF] solve highdimension-based e 1 -norm minimization problem. KCRC [START_REF] Wang | Kernel collaborative representation-based classifier for face recognition[END_REF] solves high-dimension-based e 2 -norm minimization problem. QSRC [START_REF] Zou | Quaternion collaborative and sparse representation with application to color face recognition[END_REF] solves quaternion-based e 1 -norm minimization problem. QCRC [START_REF] Zou | Quaternion collaborative and sparse representation with application to color face recognition[END_REF] solves quaternion-based e 2 -norm minimization problem. The proposed QAR solves quaternion-based e p -norm (1 < p < 2) minimization problem while the proposed HD-QAR solves high-dimension-quaternion-based e pnorm (1 < p < 2) minimization problem. [START_REF] Zou | Quaternion collaborative and sparse representation with application to color face recognition[END_REF] also compares difference between the quaternion-based sparse methods with the group-based sparse representation classification (GSRC) methods [START_REF] Huang | The benefit of group sparsity[END_REF]- [START_REF] Lv | The group lasso for stable recovery of blocksparse signal representations[END_REF], joint-based sparse representation classification (JSRC) methods [START_REF] Lee | Compressive diffuse optical tomography: noniterative exact reconstruction using joint sparsity[END_REF], [START_REF] Kim | Compressive music: Revisiting the link between compressive sensing and array signal processing[END_REF].

V. EXPERIMENT RESULTS

This section assesses the effectiveness of the proposed methods. We compare the proposed methods with CRC, ECRC-M, ECRC-A, QCRC, KCRC, SRC, ESRC-M, ESRC-A, GSRC, JSRC, QSRC and KSRC on several well-known color face databases.

A. Face Recognition with pose, lighting and expressions variations

AR face database includes 126 people (70 men and 56 women) with more than 4,000 color face images. These face images have different illumination conditions, expressions and occlusions (sun glasses and scarf). Each person has two sessions (2 different days). The same face images were taken in both sessions. The face images were taken at the CVC under strictly controlled conditions. Participants may have different wear (clothes, glasses, etc.), make-up, hair style, etc. Following [START_REF] Zou | Quaternion collaborative and sparse representation with application to color face recognition[END_REF], we select a subset including the images with illumination and expression variations for this experiment. For each person, we randomly choose n (1, 2, 3, 4, 5) face images from session 1 as training set and the face images in session 2 are utilized as testing set. To obtain the convincing results, we randomly select the images and run this process ten times. The average recognition rates of this experiment is described in Fig. 4 (a), and the mean recognition rates of all n (1, 2, 3, 4, 5) are listed in Table IV. Caltech Faces Database has 27 subjects with 450 color face images under various illumination, expression and backgrounds. Following [START_REF] Zou | Quaternion collaborative and sparse representation with application to color face recognition[END_REF], we use the Viola-Jones face detector [START_REF] Viola | Robust real-time face detection[END_REF] to detect faces. Then, n (1, 2, 3, 4, 5) face images of each person are selected as training set and the rest face images are used as testing set. To get the convincing results, we also randomly select the images and run this process ten times. Fig. 4 (b) shows the average recognition rates of this experiment. Table V lists mean recognition rates of all n (1, 2, 3, 4, 5). Color FERET (Facial Recognition Technology) Database has 119 subjects with 14, 126 face images. Under a semicontrolled environment, these images are collected in 15 sessions. Following [START_REF] Zou | Quaternion collaborative and sparse representation with application to color face recognition[END_REF], the frontal face images of the sets with the letter code fa and fb are utilized for experiments. The Viola-Jones face detector [START_REF] Viola | Robust real-time face detection[END_REF] is used to detect faces. In the experiment, we select p (20, 30, 40, 50, 60) percent of face images of each person as training set, and the rest face images are used as testing set. To get the convincing results, we also randomly select the images and run this process ten times. The average recognition rates of this experiment are shown in Fig. 4 (c). The mean recognition rates of all p (20, 30, 40, 50, 60) are listed in Table VI. Considering the Tables IV, V and VI, we conclude three points as follows. • Due to the fact that the quaternion can obtain the correlation information among different channels of color images, the quaternion-based sparse methods obtain the better performance compared to the non-quaternion-based sparse methods. • The high dimension feature in kernel space can help the sparse-based methods obtain better performance. Thus, the kernel-based sparse methods outperform the corresponding linear-based sparse methods. • The adaptive e p -norm (1 ≤ p ≤ 2) minimization can benefit from both e 1 -norm minimization and e 2 -norm minimization according to the correction of the samples in the prototype set X. Therefore, the proposed methods obtain better performance than QCRC and QSRC in [START_REF] Zou | Quaternion collaborative and sparse representation with application to color face recognition[END_REF].

B. Face Recognition with low resolution

In this section, we use the SCface (Surveillance Cameras Face) to assess the effectiveness of the proposed methods for face recognition with low resolution. The SCface database is a dataset of static face images of human. Face images were collected in uncontrolled indoor environment with 5 video surveillance cameras of various qualities. This dataset includes 130 persons with 4160 face images (in visible and infrared spectrum) from different quality cameras mimic the real-world

C. Face Recognition with complicated background

This section uses the CMU PIE face database to evaluate the effectiveness of the proposed methods for face recognition with complicated background. CMU PIE database includes more than 40,000 facial images of 68 people. Using the CMU (Carnegie Mellon University) 3D Room, this dataset is collected by imaging each person across 13 different poses, under 43 different illumination conditions, and with four different expressions. This dataset is quite challenging because all face images have the various and complicated background. The experiment is implemented as follows: one-third images of expressions of each person are selected to form training set, while the rest two third images are considered as test images. The experiment results of comparison methods are listed in Table VIII. We can observe that the proposed methods obtain the better performance than other state-of-the-art methods.

D. Face Recognition in the Wild

This section uses the LFW (Labeled Faces in the Wild) to assess the performance of the proposed methods for face recognition against environment in the wild. LFW face database has 5749 persons with 13233 color face images taken in real and unconstrained environments. For each person, the face images have large visual variations in illumination, pose, occlusion, etc. This database is quite challenging because the face images of each person have the large variations in unconstrained environments. Following [START_REF] Zou | Quaternion collaborative and sparse representation with application to color face recognition[END_REF], we choose the persons with over n (40,60, 80) face images in this experiment. Then we randomly select 10 percent of the face images per person as training set and the rest images are used as testing set. To get the convincing results, we also randomly select the images and run this process ten times. The average recognition rates of the comparison methods are listed in Fig. 5 (a). Observing this figure, we know that the proposed methods could improve the performance of other comparison methods and obtain the higher recognition rate, in spite of the face images having the large variations in unconstrained environments.

E. Face recognition with noise image

In this subsection, we use the LFW (Labeled Faces in the Wild) face database to test the robustness of the proposed methods against noise. The training sets are same as those of the above section. The testing set is obtained by Matlab function "imnoise" to insert the Gaussian white noise. Notice that only the testing face images are inserted with the Gaussian white noise, and the training face images are utilized the original images of the database. To get the convincing results, we also randomly select the images and run this process ten times. The average recognition rates of the comparison methods are listed in Fig. 5 (b). Observing this figure, we know that the proposed methods could outperform other comparison methods and obtain the higher recognition rates for face images with the noise and large variations in unconstrained environments.

VI. CONCLUSION

In this paper, quaternion-based adaptive representation (QAR) was proposed for color face recognition. QAR can adaptively combine the e 1 -norm and e 2 -norm and benefit from both. Moreover, we proposed the high-dimension quaternionbased adaptive representation (HD-QAR), which can obtain the high dimension feature information in the kernel space, Thus, HD-QAR can obtain the correlation information among the different channels. Experimental results showed that the proposed methods achieve better performance compared to several state-of-the-art methods. The analyses and the experimental results on three databases confirmed the effectiveness of the proposed methods for color face recognition.

Fig. 1 .

 1 Fig. 1. Motivation: Because quaternion-based e 2 -norm and e 1 -norm both have their own advantages, we want to design a quaternion-based ep-norm (1 ≤ p ≤ 2) that can adaptively combine quaternion-based e 2 -norm and e 1 -norm and benefit from both. The above is the geometric interpretation of ep-norm, e 2 -norm and e 1 -norm.

Algorithm 1

 1 Proposed Quaternion-based Adaptive Representation (QAR) Method Require: The original training set X with L samples and a testing sample x. Initialize the value of Z, z, α, m 1 , M 2 , u, ρ, ε and u max Ensure: Label of the testing sample x. 1: Constitute the objective function using the original training set X and testing sample x as ℵ( α) * = arg min ||ℵ( ẋ) -( Ẋ)ℵ( α)|| 1 + λ|| ( Ẋ)Diag(ℵ( α))|| * 2: Transform the objective function to a Lagrange multiplier problem as

9 :Fig. 2 .

 92 Fig. 2. Motivation of HD-AQR. Given a color face image, it has red, green and blue channel. Quaternion transformation can obtain the correlation information among different color channels. The kernel quaternion transformation can obtain the high dimension correlation information among different color channels.

2 : 3 :

 23 Use K and k to constitute the objective function as min S,s,γ ||S|| * + ||s|| 1 s.t. s = k-Kγ S = KDiag(γ) Transform the objective function to a Lagrange multiplier problem as

10 :

 10 Check convergence by ||k -Kγ -s|| ∞ ≤ ε and ||S -KDiag(γ)|| ∞ ≤ ε 11: end while 12: Compute the distance d c = ||k-Kcαc|| ||αc|| and classify the testing sample x by c * = arg min(d c ) ed as. min || φ ( Ẋ)Diag(γ)|| * s.t. ℵ φ ( ẋ) = φ ( Ẋ)γ

Fig. 3 .

 3 Fig. 3. Several color images from the (a) AR database, (b) Caltech database, (c) colorFERET database, (d) SCface database, (e) CMU-PIE database and (f) LFW database.

  min ||α|| 1 s.t. φ(x) = φ(X)α high-dimension-based e 1 -norm minimization KCRC min ||α|| 2 s.t. φ(x) = φ(X)α high-dimension-based e 2 -norm minimization QSRC min ||α|| 1 s.t. ẋ = Ẋα quaternion-based e 1 -norm minimization QCRC min ||α|| 2 s.t. ẋ = Ẋα quaternion-based e 2 -norm minimization QAR min || Ẋα|| * s.t. ẋ = Ẋα quaternion-based ep-norm (1 < p < 2) minimization HD-QAR min ||φ( Ẋ)α|| * s.t. φ( ẋ) = φ( Ẋ)α high-dimension-quaternion-based ep-norm (1 < p < 2) minimization solved by the following minimization problem γ * = arg min γ L(S, s, γ)

Fig. 4 .

 4 Fig. 4. The recognition rates of several methods on the (a) AR database, (b) Caltech database, (c) colorFERET database and (d) SCface database.

Fig. 5 .

 5 Fig. 5. The recognition rates of several methods on the (a) LFW face database with normal image and (b) LFW face database with noise image.

  TABLE III IN[START_REF] Zou | Quaternion collaborative and sparse representation with application to color face recognition[END_REF]

	Method	QCRC	QSRC
	1	48.70	51.57
	2	67.94	68.67
	3	75.80	75.41
	4	78.99	78.41
	C. Representation-based Methods
	Wright et al. proposed sparse representation classification
	(SRC) in 2009. It represents the testing sample by the linear
	combination of training samples of all classes. SRC solves
	e 1 -minimization optimization problem as
	min ||α|| 1 s.t. x = Xα	(8)
	To reduce the computation cost and obtain the better represen-

tation, Zhang et al. proposed the collaborative representation based classification (CRC). In CRC, the authors argued that the collaborative representation should be better than the linear combination of training samples. CRC solves the e 2minimization optimization problem as min ||α|| 2 2 s.t. x = Xα

TABLE III THE

 III PROPOSE METHODS VS WELL-KNOWN SPARSE METHODS.

Methods

Objective Function Solved minimization problem SRC min ||α|| 1 s.t. x = Xα e 1 -norm minimization CRC min ||α|| 2 s.t. x = Xα e 2 -norm minimization KSRC

TABLE IV THE

 IV RECOGNITION RATES OF SEVERAL METHODS ON AR FACE DATABASE

	Methods	Mean Recognition Rate
	CRC	63.25
	KCRC	63.81
	ECRC-M	63.81
	ECRC-A	68.74
	QCRC	70.53
	SRC	58.21
	KSRC	58.84
	ESRC-M	66.95
	ESRC-A	67.65
	GSRC	69.89
	JSRC	69.56
	QSRC	71.03
	QAR	72.42
	HD-QAR	72.87
		TABLE V
	THE RECOGNITION RATES OF SEVERAL METHODS ON CALTECH FACE
		DATABASE
	Methods	Mean Recognition Rate
	CRC	83.60
	KCRC	84.09
	ECRC-M	85.40
	ECRC-A	86.00
	QCRC	88.87
	SRC	83.07
	KSRC	83.47
	ESRC-M	86.77
	ESRC-A	86.20
	GSRC	86.68
	JSRC	86.95
	QSRC	88.52
	QAR	90.29
	HD-QAR	90.70

TABLE VI THE

 VI RECOGNITION RATES OF SEVERAL METHODS ON COLORFERET FACE

		DATABASE
	Methods	Mean Recognition Rate
	CRC	70.24
	KCRC	70.76
	ECRC-M	72.17
	ECRC-A	72.49
	QCRC	74.82
	SRC	70.73
	KSRC	71.27
	ESRC-M	72.92
	ESRC-A	74.01
	GSRC	73.77
	JSRC	73.43
	QSRC	75.14
	QAR	76.97
	HD-QAR	77.39
	conditions. The sizes of face images are various. Due to some
	of the surveillance face images are of quite low resolution
	and quality, this dataset is quite challenging. Following [17],

we randomly choose 2080 face images of 130 persons in this experiment. Then, p (10, 20, 30, 40, 50) percent of face images of each person are used as training set, and the rest face images are used as testing set. To achieve the convincing results, we also randomly select the images and run this process ten times. The average recognition rates of this experiment are plotted in Fig.

4 (d

) and the mean recognition rates of all p (10, 20,

TABLE VII THE

 VII RECOGNITION RATES OF SEVERAL METHODS ON SCFACE FACE

		DATABASE
	Methods	Mean Recognition Rate
	CRC	62.46
	KCRC	63.16
	ECRC-M	63.18
	ECRC-A	64.47
	QCRC	65.53
	SRC	60.96
	KSRC	61.32
	ESRC-M	62.02
	ESRC-A	62.48
	GSRC	63.50
	JSRC	63.53
	QSRC	64.16
	QAR	67.15
	HD-QAR	67.72
	30, 40, 50) are listed in Table VII. From this table, we can get
	the following conclusions. Quaternion-based sparse methods

outperform the non-quaternion-based sparse methods. Kernelbased sparse methods have better performance than linearbased sparse methods. The adaptive e p -norm (1 < p < 2) minimization can benefit from both e 1 -norm minimization and e 2 -norm minimization according to the correction of the samples in the prototype set X. Therefore, the proposed methods obtain the better performance than QCRC and QSRC

TABLE VIII THE

 VIII RECOGNITION RATES OF SEVERAL METHODS ON CMUPIE FACE

		DATABASE
	Methods	Mean Recognition Rate
	CRC	75.12
	KCRC	76.10
	ECRC-M	79.42
	ECRC-A	78.10
	QCRC	83.84
	SRC	77.27
	KSRC	78.04
	ESRC-M	84.12
	ESRC-A	83.79
	QSRC	85.04
	QAR	86.27
	HD-QAR	86.79

https://en.wikipedia.org/wiki/Matrix norm

This work was supported in part by the Macau Science and Technology Development Fund under Grant FDCT/016/2015/A1 and by the Research Committee at University of Macau under Grants MYRG2014-00003-FST and MYRG2016-00123-FST. (Corresponding author is Yicong Zhou.)


Lesage Sylvain
email: sylvain.lesage@irisa.fr

Gribonval Rémi
email: remi.gribonval@irisa.fr

Bimbot Frédéric
email: frederic.bimbot@irisa.fr

Benaroya Laurent

PREPRINT -LEARNING UNIONS OF ORTHONORMAL BASES WITH THRESHOLDED SINGULAR VALUE DECOMPOSITION

We propose a new method to learn overcomplete dictionaries for sparse coding. The method is designed to learn dictionaries structured as unions of orthonormal bases. The interest of such a structure is manifold. Indeed, it seems that many signals or images can be modeled as the superimposition of several layers with sparse decompositions in as many bases. Moreover, in such dictionaries, the efficient Block Coordinate Relaxation (BCR) algorithm can be used to compute sparse decompositions. We show that it is possible to design an iterative learning algorithm that produces a dictionary with the required structure. Each step is based on the coefficients estimation, using a variant of BCR, followed by the update of one chosen basis, using Singular Value Decomposition. We assess experimentally how well the learning algorithm recovers dictionaries that may or may not have the required structure, and to what extent the noise level is a disturbing factor.

INTRODUCTION

Sparse coding [START_REF] Lewicki | Learning overcomplete representations[END_REF][START_REF] Lewicki | A probabilistic framework for the adaptation and comparison of image codes[END_REF] is a useful tool to analyze and try to explain the structure of series of observed data, such as successive time frames of an audio signal [START_REF] Abdallah | If edges are the independent components of natural images, what are the independent components of natural sounds?[END_REF] or natural images [START_REF] Bell | The 'independent components' of natural scenes are edge filters[END_REF]. Formally, assume that we observe T vectors x(t) = (x n (t)) N n=1 , 1 ≤ t ≤ T which are supposed to have been generated following the model:

x(t) = As(t) + ϵ(t) (1)
A being an overcomplete dictionary (an N ×K matrix, with K ≥ N ), s(t) ∈ R K some "sparse" coefficients and ϵ(t) ∈ R N a Gaussian noise. Sparse coding can be viewed as a way of estimating A from the only observation of X = AS + E where X is the N × T matrix containing T signal frames (similar notations being used for S and E).

Jointly optimizing the coefficients and the dictionary, under constraints added to enforce the well-posedness of the problem, is a hard task, so we use an alternating optimization strategy:

⋆ The initial part of this work was done in collaboration with Laurent Benaroya while he was finishing his PhD with IRISA. [START_REF] Lewicki | Learning overcomplete representations[END_REF]. Coefficient update given a dictionary A:

arg min S ∥X -AS∥ 2 2 + λ∥S∥ 1 (2)
2. Dictionary update given coefficients S:

arg min A ∥X -AS∥ 2 2 (3)
under some constraint on A.

The coefficient update step (2) can be justified in a probabilistic framework using a Laplacian prior on the coefficients s k (t) [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]. Moreover, it is a simpler parent problem of the NP-hard combinatorial problem, where ∥S∥ 1 is replaced with ∥S∥ 0 0 , the number of non-zero components in S. Computing the solution to Eq. ( 2) by Quadratic Programming is rather computationally intensive in the general case where A has no special structure. However, when A is a union of orthonormal bases (ONB), Block Coordinate Relaxation (BCR) methods are efficient [START_REF] Sardy | Block coordinate relaxation methods for nonparametric signal denoising with wavelet dictionaries[END_REF]. Another motivation to constrain the dictionary to be a union of ONB is that it seems that audio signals [START_REF] Molla | Determining local transientness in audio signals[END_REF] and images [START_REF] Starck | Image decomposition via the combination of sparse representations and a variational approach[END_REF] can indeed be modeled as the superimposition of several layers, each of which having a sparse representation in its own adapted ONB. Note that when A is constrained to have this precise structure, the dictionary update step (3) is also made relatively easy. This step, in a probabilistic framework, can be interpreted as a likelihood maximization and solved with an Expectation-Maximization (EM) algorithm [START_REF] Moon | The expectation-maximization algorithm[END_REF].

In Section 2, we describe BCR and a variant which we used to solve [START_REF] Lewicki | A probabilistic framework for the adaptation and comparison of image codes[END_REF]. In Section 3, we introduce our algorithm to learn a union of bases by iteratively optimizing (3) with respect to each basis. In Section 4 we describe and analyze the experiments made with the learning algorithm on data generated following the model [START_REF] Lewicki | Learning overcomplete representations[END_REF]. We study the influence of the number T of frames of the learning dataset, of the a priori knowledge on the noise level, and of the possible modeling error corresponding to the fact that the true A might not be a union of bases or the number of bases could be wrong.

COMPUTATION OF SPARSE COEFFICIENTS

Finding sparse coefficients for the observed data X is the result of a compromise between

• the minimization of the reconstruction error:

∥X -AS∥ 2 2 := N n=1 T t=1 |x n (t) -(As(t)) n | 2 .
• the minimization of a diversity measure, the most common ones being:

∥S∥ p p := K k=1 T t=1 |s k (t)| p for 0 ≤ p ≤ 1.
The strict diversity, defined by the number of non-zeros coefficients, is given by ∥S∥ 0 0 .

This problem is generally difficult. It is indeed NP-hard with the ∥S∥ 0 0 diversity measure when A is an arbitrary redundant dictionary. Many sub-optimal algorithms have been proposed such as Matching Pursuit (MP) [START_REF] Mallat | Matching pursuits with timefrequency dictionaries[END_REF], Basis Pursuit (BP) [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] and FOCUSS [START_REF] Kreutz-Delgado | Dictionary learning algorithms for sparse representation[END_REF]. The latter solves the minimization problem

min S ∥X -AS∥ 2 2 + λ∥S∥ p p (4)
and Basis Pursuit solves it for p = 1 (see Eq. ( 2)). These algorithms are generally rather computationally intensive. However Basis Pursuit can be implemented more efficiently with a Block Coordinate Relaxation (BCR) method [START_REF] Sardy | Block coordinate relaxation methods for nonparametric signal denoising with wavelet dictionaries[END_REF] when A is a union of ONB. Moreover, it has been shown that, under some conditions, the solution of (4) for any 0 ≤ p ≤ 1 is close to the solution given by Basis Pursuit [START_REF] Gribonval | On the strong uniqueness of highly sparse representations from redundant dictionaries[END_REF].

In this section we recall how Basis Pursuit is implemented with soft-thresholding when A is a single ONB, then we remind the reader about BCR and eventually we describe a variant of BCR that we introduced to deal with low noise levels (small threshold parameter λ).

Case of an orthonormal basis

When A is an orthonormal basis, the solution of ( 2) is given by soft thresholding:

∀k, t ŝk (t) = ⎧ ⎨ ⎩ a T k x(t) -λ/2 if a T k x(t) > λ/2 0 if |a T k x(t)| ≤ λ/2 a T k x(t) + λ/2 if a T k x(t) < -λ/2
(5) where a k is the k th column of A (also called atom of the dictionary).

Case of a union of orthonormal bases

When A = [A 1 , . . . , A L ] is a union of L orthonormal bases, the coefficients S are decomposed in L subsets of coefficients S l corresponding to the L bases, as S = [S T 1 , . . . , S T L ] T . The BCR algorithm, described in Table 1 deals with the difficulty to directly solve (2) for a redundant dictionary A by successively solving it for its different bases A l . Then the sub-coefficients of an initial estimate S init are iteratively updated until convergence is reached. The BCR algorithm has been proven to converge in a wake sense to a solution of (2), when, in Step 1, the selection of S l follows a systematic cycle, or results from an optimal descent rule [START_REF] Sardy | Block coordinate relaxation methods for nonparametric signal denoising with wavelet dictionaries[END_REF]. Unfortunately, if λ is very small (which corresponds to looking for a small reconstruction error, namely the low noise assumption), BCR might converge very slowly since almost no thresholding is performed in Step 3. In order to compute sparse coefficients in this low noise case, we propose to start BCR with a large initial threshold λ 0 and decrease it regularly, leading to the algorithm explained in Table 2. During the very first steps, since the threshold is high, sparsity is enforced; when the threshold becomes smaller, the error vanishes.

1. Select a subset S l of the current S to update;

2. Compute X l = X -i̸ =l A i S i ;
3. Update S l by replacing it by arg min

S l ∥X l -A l S l ∥ 2 2 + λ∥S l ∥ 1 ,
which is computed by soft thresholding (Eq. ( 5));

4. If the stopping criterion is not reached, go to step 1.

Experiments

Even though we have no proof of convergence of this modified BCR algorithm, we observed experimentally that if the two parameters N it and λ 0 are appropriately chosen, it reaches a solution close to those provided at a higher computational cost by MP and FOCUSS. We ran the following experiment. Data are sparsely generated from a dictionary that is the union of two orthonormal bases. Using this dictionary, the BCR variant was run on the data, for different values of the initial threshold (between 10 -6 and 10 3 ), and for different number of iterations (3, 10 and 100). The diversity of the obtained coefficients is displayed in figure 1. We can see that, for any N it , the initial threshold λ 0 greatly impacts the diversity of the estimated coefficients. An optimal value of λ 0 may be choosen a posteriori to minimize [START_REF] Lewicki | A probabilistic framework for the adaptation and comparison of image codes[END_REF]. Note that the higher the number of iterations, the sparser the coefficients obtained using the optimal threshold.

For N it = 100 and λ 0 in a fairly large range, the diversity of the coefficients obtained by the BCR variant is the same than for MP and FOCUSS. However, the BCR algorithm variant is computationally less costly, taking about 3.5 seconds instead of about 9 seconds for FOCUSS, and 10 seconds for MP.

DICTIONARY LEARNING WITH SVD

The algorithm used to learn a union of L orthonormal bases (ONB) is described in Table 3. To understand the rationale behind the use of the SVD in Step 3, we will first analyze the optimization problem (3) when A is constrained to be a single ONB. Then, we will briefly explain how the algorithm for L bases is derived from the single basis one, and we will discuss in more details how the coefficient update (Step 2) is performed, depending whether we know which value of λ to use in (Step 2) or we want to adapt it to the data. As for the the stopping criterion (Step 4), we simply set a priori the number of learning steps. Studying how much the dictionary varies between two steps may help design a better criterion in the future.

Choose an initial dictionary

A = [A 1 . . . A L ];
2. Update the coefficients S = [S 1 . . . S L ] using the current A (see text);

3. Choose which basis A l to update and:

(a) Compute X l = X -i̸ =l A i S i
(b) Compute a singular value decomposition:

X l S T l = UDV T (c) Update A l = UV T 4.
If the stopping criterion is not reached, go to step 2 (see text).

Table 3. Learning algorithm for L orthonormal bases

Learning a single orthonormal basis

The optimization problem (3) with the constraint that A is an ONB can be written as the minimization of a Lagrangian,

L(A, µ) = ∥X -AS∥ 2 2 + Tr µ(A T A -Id) (6)
where µ is an N ×N matrix of Lagrange multipliers chosen so that the minimizing matrix A opt is an ONB. The optimal dictionary A opt , minimizing the Lagrangian, is obtained as follows (the proof is in annex A): Proposition 3.1 Let UDV T be (one of) the Singular Value Decomposition(s) (SVD) of XS T , that is to say U and V are orthogonal matrices and D is a diagonal matrix. An optimal ONB is given by A opt = UV T . If XS T is invertible, it is the unique optimal solution.

Note that, when XS T is invertible, the product UV T does not depend on the particular choice of a SVD UDV T .

Learning a union of L orthonormal bases

Ideally, when A is constrained to be the union of L orthonormal bases, one would like to perform the dictionary update step (3) by minimizing the Lagrangian:

∥X - L l=1 A l S l ∥ 2 2 + L l=1 Tr µ l (A T l A l -Id) .
However, this optimization problem does not have an explicit solution as in the case of a single ONB. The principle behind the algorithm described in Table 3 is that, at each iteration, only one of the bases A l is optimized.

If we know which parameter λ to use in Eq. ( 2) -for example, if we know the prior distributions of ϵ and s in the probabilistic model (1) -then we perform Step 2 with the regular BCR algorithm. In many practical cases however, it is difficult to have an idea of a relevant value for λ, or in the low noise limit, λ is too small. Below we discuss how Step 2 in Table 3 is performed depending on the noise model.

Unknown gaussian noise

If Gaussian noise ϵ(t) is assumed, with unknown variance, we use the algorithm proposed by Azzalini et al, [START_REF] Azzalini | A recursive algorithm for nonlinear wavelet thresholding : Application to signal and image processing[END_REF]: starting from an exact representation X = AC, with A the current dictionary, and the estimate S = 0 we iterate the following steps 1. compute the variance σ 2 of the residual R = C -S; 2. update S by letting it contain all the entries of R that are above the threshold λ = 2 log(N )σ 2 .

3. if the last update did not modify S, then stop, else go to 1.

To compute the exact decomposition C, rather than using C = A + X, one can use indifferently FOCUSS, MP, or the variant of BCR, because they all encourage sparsity of the coefficients.

Low noise limit

In the low noise limit, each of the above strategies for Step 2 will essentially fail since the dictionary update step will (almost) not change the dictionary. To see it, let us simply analyse the first iteration of the algorithm described in Table 3.

Given the initial dictionary A = [A 1 . . . A L ] and the data X, the coefficients S = [S 1 . . . S L ] = S(A, X) computed at Step 2 are uniquely defined. Because of the low noise assumption, they give a perfet reconstruction X = AS.

Then, if X l S T l is invertible, the proposition 3.1 shows that the optimal basis Âl which minimizes the Lagrangian (6) is unique. And, for the current basis A l , the first term in the Lagrangian -the error term -is null, and the second term is null too, because of the orthonormality of the basis. The lower bound of the Lagrangian for an orthonormal basis is then reached by the current basis. We conclude that Âl = A l and there can be no change.

If XS T is not invertible, the optimal basis may be different, but this case does not appear experimentally.

Thus, in the low noise limit, we propose to perform Step 2 with a heuristic strategy that allows the dictionary to be modified at each step. To do that, we add some reconstruction error at each step by thresholding the coefficients. The added reconstruction error is decreased each time Step 2 is performed so that at the end of the process we get exact reconstruction.

1. Choose an initial dictionary A = [A 1 . . . A L ] and set N kept , the number of kept coefficients by frame, to zero;

2. Update the coefficients S = [S 1 . . . S L ] using the current A, by the BCR variant;

3. If the algorithm stagnates, increment N kept ;

4. Set all but the N kept T greatest coefficients of S to zero, giving S thres ;

5. Choose which basis A l to update and:

(a) Compute X l = X -i̸ =l A i S thres i (b)
Compute a singular value decomposition:

X l S thres T l = UDV T (c) Update A l = UV T 6.
If the stopping criterion is not reached, go to step 2.

Table 4.

Learning algorithm for L orthonormal bases in the low noise limit Since the learning (experimentally) leads to decrease the diversity of the coefficients S along the iterations, the algorithm is said to stagnate when this diversity stops decreasing. By example, it can be decided the algorithm stagnates when the relative difference between the diversities of the five last coefficients matrices is least than 10 -3 .

EXPERIMENTS

The following experiments have been designed with synthetic data generated with the model (1) using a known reference dictionary A ref . The goal is to study the influence of various parameters on the performance, and to see how modeling error could also impact the results. We use two relevant and complementary performance measures: the false alarm rate and the missed detection rate, corresponding respectively to the relative number of estimated atoms a est (i.e. columns of A est , the dictionary estimated using the learning algorithm) that "do not match" any reference atom a ref , and to the relative number of reference atoms that "are not matched" by any estimated atom. In order to evaluate the performance on a wide range of conditions, each experiment is run with N r different dictionaries, and the performance measures are averaged over these runs.

Influence of the number T of signal frames

First, we study the impact of the number T of frames used to learn the dictionary, for different dimensions N . Data are generated with A ref a union of two random ONB, using the noiseless model (1) with s(t) (of dimension 2N ) containing between 1 and r, randomly located, non-zero coefficients, that follow a standard Gaussian law. r is depending on the dimension N and is chosen as follows: r = ⌊(1 + √ N )/2⌋. Indeed, for this number of non-zero coefficients, the generating s(t) are the only coefficients, among those giving exact decomposition, that solve (2) [START_REF] Gribonval | On the strong uniqueness of highly sparse representations from redundant dictionaries[END_REF].

The experiments are performed for the following dimensions : N = 4, 8, 16, 32, 64, and, for each of them, for T = N/4, N, 2N, 5N, 10N, 20N, 50N, 100N . The MD Rate is decreasing as the relative number of learning frames is growing. For T ≤ 2N , the dictionary estimation yields poor performance, the false alarm rate and the missed detection rate are both greater than 50%, for ξ = 0.99, and even than 90%, for N ≥ 16. On the contrary, for T = 100N , about seven out of ten atoms are correctly estimated. Noting that the computing time greatly increases with T , we set T = 50N in the rest of the experiments.

Influence of the noise level

To understand the effect of the noise level on the performance, we repeat the above experiments with the same data to which we add noise at various signal to noise ratios (SNR): +∞ dB, 10 dB, 0 dB. At each SNR level we run the three configurations of the learning algorithm designed respectively for known λ, small λ and unknown λ.

Learning algorithm +∞ dB +10 dB +0 dB λ known 7% 28% 59% λ small 6% 30% 63% λ unknown 42% 58% 86%
Table 5. Missed Detection Rate depending on the noise level on data, and on the learning method, with ξ = 0.99 (average over N r = 60 experiments) Table 4.2 shows that the algorithms with known λ, and small λ (without prior knowledge on λ), give almost similar results, the first one performing better when there is some noise, while the second one giving better estimation in the low-noise case. Unfortunately, the algorithm designed for unknown λ always misses more atoms. Thus, in the following, this algorithm by Azzalini will no more be used.

The three configurations of the algorithm are greatly dependent on the noise level. Indeed, even for a +10dB signal to noise ratio, among the N r = 60 runs, the Missed Detection Rate of the most succesful experiment is no smaller than 14%. This means that for this SNR level, the algorithms never finds more than 86% of the atoms of the dictionary.

On the contrary, when there is no noise, the first two algorithms retrieve all the atoms in most of the runs. A part of the experiments entirely fail, leading to the means displayed in Table 4.2 (7% and 6%).

Influence of the dimension N of signal frames

In order to understand the influence of the dimension N of signal frames on the behaviour of the algorithm, we run the learning for noiseless data, as in subsection 4.1, for different values of N . The dictionary is made of two random ONB.

The first experiment illustrates how many iterations are required for the algorithm to converge, depending on N . For each N , N r = 100 experiments are launched, and results are then averaged. The algorithm is decided to converge when the original dictionary has been retrieved with a sufficient precision, namely when the Missed Detection Rate, for a threshold ξ = 0.99, reaches 0%. The size of frames N is chosen between 4 and 128. Indeed, for higher values of N , the computation is too long, and learning is unfeasible.

The figure 3 shows the number of iterations needed to converge, depending on N , when the algorithm converges. The N r points, for each N , represent the number of iterations needed to converge in each of the N r experiments. The circles represent the mean for a given N , and the bars, the associated standard deviation. The broken line is the linear regression of the means. We see that, in the case when the learning is succesful, the average number of iterations to reach convergence increases linearly with the size N . Using the mean is accurate since the variance of the data points around the mean is small enough to clearly see the linear dependence between N and the number of iterations before convergence.

The figure 4 shows the percentage of converging experiments.

We see that, for very small values of N , very few experiments converge, and the algorithm seems also to be less efficient when N become high. Work should be done to decide, for a given experiment, if the algorithm has converged to a correct solution. An idea would be to learn three dictionaries and to make a vote, hoping that two dictionaries out of the three are the good ones.

A second experiment is performed, showing the influence of the size of the signal frames on the computational cost. The number of samples N total in the whole signal is set constant (i.e. for different N , the signals are different frame hashings of the same original signal). The number T of frames is then T = N total N . Note we could expect bad learning results when T become too small, according to Percentage of experiments where the algorithm have converged (the estimated dictionary is the original one) subsection 4. In fact, in this experiment, we only look at the computational time, and learning is not performed until the end (and the result of the learning is not important here).

The average computational time of each step of the algorithm in Table 3 is measured. As a result, we observe that for high values of N (from N = 128 to 2048), the computational time of the SVD Step 3b, the matrices multiplication Steps 3a and 3c and the coefficients update Step 2 grow as a power of N (between O(N 2 ) and O(N 3 )). Thus, on a Pentium 4, for N = 512, a single iteration takes 10s, and for N = 2048, it takes 324s.

The computational time for performing learning then increases for two reasons as N grows. Firstly, it increases due to the higher computational cost of each step, and secondly, more iterations are needed before convergence. The algorithm is then untractable for values of N greater than 64.

But, in audio signal, the frame size is commonly about 20ms. At a sampling frequency of 44100Hz, a frame of 23.2ms corresponds to N = 1024, that is untractable with this algorithm. A challenge would be to adapt the method to high dimension frames.

Influence of the model mismatch

We designed experiments to analyse the behavior of the learning algorithm when there is a mismatch between the number of ONB in the reference dictionary and in the estimated one.

We run the same algorithm to estimate a pair of ONB on three datasets generated as above with the noiseless model (1) with

• A ref,1 a single random ONB • A ref,2 a

union of two ONB

• A ref,3 a union of three random ONB.

Note that with A ref,3 , if the three bases are sufficiently different one from another, one cannot expect to get less than 33% missed detection, because only 2N atoms are estimated while there are 3N reference atoms. With A ref,1 , as soon as there is no more than 50% false alarm we are sure to have recovered the atoms of the reference dictionary, but they may be split in the two learned bases.

Reference dictionary

A ref,1 A ref,2 A ref,3
Average missed detection rate 0.5% 7% 99.5% Average false alarm rate 44% 7% 99% The results of two hundred runs are summarized in Table 4.4. More precisely:

• A ref,1 : almost all reference atoms are retrieved. In 55 cases out of 100, one of the estimated basis is exactly A est,1 . In the 45 other cases out of 100, the retrieved atoms are shared between the two bases, 82% in the first basis, and 17% in the second, while 1% are not detected.

• A ref,2 : the average performance values hide two distinct behaviours. In 92% of the experiments, the dictionary is perfectly estimated, while in 8% of the cases, learning totally failed, without any well estimated atom at all.

• A ref,3 : never more than 15% of the atoms were retrieved, the average being only 0.5%.

The algorithm seems to be efficient only when estimating at least as many bases as there are in the reference dictionary. A good strategy could therefore be to learn a lot of bases, and to estimate a posteriori the number of interesting ones.

CONCLUSION

We have presented a new method for learning, from a set of observed data vectors, a dictionary structured as a union of orthonormal bases, with the objective that the decomposition of the data on this dictionary would be sparse. We have demonstrated on synthetically generated data that this method is able to recover a relevant underlying dictionary provided that it knows a priori the structure (i.e. the number of ONB in the dictionary). The approach seems to behave reasonably well even when the number of bases is overestimated. We are now considering several remaining practical problems, namely estimating the number of bases, studying how the algorithm scales when the dimension N of the data becomes large and extending experiments to real audio signals or images. Last but not least, we are looking for conditions where we can prove that the algorithm converges to the underlying dictionary.

APPENDIX A. PROOF OF PROPOSITION 3.1

Let detail the proof giving the optimal solution A opt to the minimization of the Lagrangian (6). The Lagragian is:

L(A, µ) = ∥X -AS∥ 2 2 + Tr µ(A T A -Id)
The optimal dictionary A opt necessarily verify that the gradients ∇ µ L and ∇ A L are equal to zero:

A T A -Id = 0 (7) -2(X -AS)S T + A(µ + µ T ) = 0 (8) or A T A = Id XS T = A SS T + 1 2 (µ + µ T ) (9)
Let Z := XS T and Y := SS T + (µ + µ T )/2. While Z can be explicitly computed, Y is unknown since it depends on the unknown multipliers µ.

The system (9) implies:

ZZ T = AYY T A T = AY 2 A T
because Y is symmetric, and

Z T Z = Y T A T AY = Y 2
because of the orthogonality of A. Then, combining the two preceding results, A opt must satisfy:

AZ T ZA T = ZZ T (10)
Let Z = UDV T be one of the possible Singular Value Decompositions (SVD) of Z, that is to say U and V are orthogonal matrices and D is a diagonal matrix. We assume that the diagonal elements of D (the singular values) are ordered decreasingly, i.e. noting

• L the number of different singular values in D,

• N k the multiplicity of the lth singular value λ k ,

• I k the identity matrix of size N l we have λ 1 > λ 2 > . . . > λ L ≥ 0, and

D = ⎡ ⎢ ⎣ λ 1 I 1 0 0 0 . . . 0 0 0 λ L I L ⎤ ⎥ ⎦
Note that between two different SVD, only the orthogonal matrices differs, while the diagonal matrix D is the same. Then equation [START_REF] Mallat | Matching pursuits with timefrequency dictionaries[END_REF] gives:

AVD 2 V T A T = UD 2 U T
and noting R = U T AV (equivalent to A = URV T ), R is orthonormal and :

RD 2 = D 2 R (11)
Let note:

R = [R ij ] i,j
where R ij is a N i by N j submatrix.

As R commutes with D 2 in (11), we have:

RD 2 = λ 2 j R ij ij = D 2 R = λ 2 i R ij ij
and then, for all i ̸ = j, One must then maximize each of the terms of the sum. For each l such that λ l > 0, one must maximize Tr(R ll ).

λ 2 j R ij = λ 2 i R ij Since λ 2 i ̸ = λ 2 j , R ij = 0 for all i ̸ = j,
As all the diagonal elements of an orthonormal matrix are lower than 1, the maximal trace is N l . It is only reached by the identity matrix, which is the only orthogonal matrix all diagonal elements of which are ones.. The trace of R ll is then maximal if and only if R ll is the identity matrix. If the last singular value λ L is null, any orthogonal matrix R LL gives λ L Tr(R LL ) = 0, and then is suitable.

The best matrix R is then the identity matrix, and it's the only one, in the case where all the singular values are strictly positive. When the last singular values are null, the identity matrix is again optimal, but no longer unique.

The optimal dictionary A, for a given SVD is then the following

A opt = UV T (13)
Remarking that the product UV T is independent from the chosen SVD, this matrix is optimal, for every chosen decomposition.

Let us prove that UV T is independent of the SVD. JE VAIS REGARDER A LA BIBLIOTHEQUE POUR UN BOUQUIN ETABLISSANT CE RESULTAT Let UDV T and LDM T be two different singular value decompositions of the matrix Z. Then,

Z T Z = VD 2 V T = MD 2 M T
Noting T = (V T M), this is equivalent to:

TD 2 = D 2 T
As before, the matrix T is then block-diagonal. Moreover, as T is orthogonal, every diagonal block is orthogonal.

The first result is that two different SVD are related by: M = VT and L = UT with T a block-diagonal matrix of orthogonal matrices, each of size of the multiplicity of the corresponding singular value. The second result is, for the same two SVD, that the products UV T and LM T are equal:

LM T = UTT T V T = UV T
since T is an orthogonal matrix.

10

 10

Fig. 1 .

 1 Fig. 1. Diversity of the coefficients computed by the BCR variant, depending on the initial threshold, and the number of iterations

Fig. 2 .Figure 2

 22 Fig. 2. Average Missed Detection Rate (on N r = 100 runs) depending on the number of frames T for different dimensions N Figure 2 displays the average Missed Detection Rate of the learned dictionary (average on N r = 100 runs) depending on the number of frames T . The False Alarm Rate is not displayed since it is nearly the same.The MD Rate is decreasing as the relative number of learning frames is growing. For T ≤ 2N , the dictionary estimation yields poor performance, the false alarm rate and the missed detection rate are both greater than 50%, for ξ = 0.99, and even than 90%, for N ≥ 16. On the contrary,

Fig. 3 .

 3 Fig. 3. Average number of iterations before convergence, depending on N , in the cases when the algorithm effectively converge

  Fig. 4. Percentage of experiments where the algorithm have converged (the estimated dictionary is the original one)

2 =

 2 l, the block R ll is an orthogonal matrix. The functional (6) must then have the following form :L(A) = ∥X-AS∥ 2 Tr(X T X)+Tr(S T S)-2Tr(XS T A T )The first two terms are independent of A, then one must maximize Tr(XS T A T ). Rewriting it givesTr(XS T A T ) = Tr(ZA T ) (12) = Tr(UDR T U T ) = Tr(DR T ) = Tr(RD) = L l=1 λ l Tr(R ll )

Table 1 .

 1 Block

	Coordinate Relaxation algorithm
	for it = 0 to N it
	use BCR with threshold λ 0 (1 -it Nit ) to update S
	end

Table 2 .

 2 Modified BCR algorithm for the low noise case

  Since all atoms have unit norm, a est and a ref are considered to match if their inner product |a T est a ref | is close enough to one. So we use a parameter ξ to decide that they match if and only if |a T est a ref | ≥ ξ. Different values of ξ yield different but related performance measures.

Table 6 .

 6 missed detection rate and false alarm rate (ξ = 0.99) depending on the number of bases in A ref , when the estimated dictionary owes two bases (average over N r = 200 runs)